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ABSTRACT 


This  report  describes  the  application  of  hybrid  computing  techniques  to 
three  areas  of  scientific  investigation.  The  areas  of  investigation 
were : 

I.  A  Hybrid  Computer  Air  Density  Model  for  a  Satellite 
Trajectory  Program 

II.  Propagation  of  Electromagnetic  Energy  in  the  Ionosphere 

III.  An  Optically  Scale’  Nuclear  Emulsive  Track  Tracer 

The  text  comprises  three  separate  sections,  or  reports.  Each  report 
includes  a  discussion  of:  mathematical  modeling,  considerations  for 
selecting  the  hybrid  approach,  hybrid  computer  implementation,  conclusions 
and  recommendations. 
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A  HYBRID  COMPUTER  AIR  DENSITY 
MODEL  IMPLEMENTATION  FOR  A 
SATELLITE  TRAJECTORY  PROGRAM 


1.1.0  INTRODUCTION 


The  ,':.rpose  of  this  work  was  to  modify  the  previously  existing 
orbital  satellite  simulation  (Contract  No.  F19628-67-C-0359)  to 
include  a  more  sophisticated  air  density  model  implemented  on  the 
hybrid  computer.  Originally,  air  density  was  provided  by  a  digital 
subroutine  which  modeled  it  as  being  a  function  only  of  altitude. 

This  over  simplification  is  unrealistic  in  that  it  neglects  many 
other  variables  that  have  a  signficant  effect  on  the  density. 

Among  the  more  obvious  of  these  are  latitude  and  longitude,  and  the 
relative  position  of  the  sun. 

The  new  density  model  that  was  implemented  is  taken  from  U.  S. 
Standard  Atmosphere  Supplements,  1966,  Part  3.  Latitude,  longitude, 
hour  angle  of  the  sun,  da}'  of  the  year,  season,  geomagnetic  activity, 
and  solar  activity  effects  are  included  in  several  formulas  which 

determine  the  exospheric  temperature,  ,  i.e.  the  temoerature 
at  the  outermost  portion  of  the  atmosphere. 

T  is  used  as  a  term  in  solving  a  set  of  differential  equations 

a>  3 

which  yield  the  number  density,  n, ,  (molecules  per  cm  )  for  each 

significant  atmospheric  constituent,  as  a  function  of  altitude. 

The  n^,  along  with  their  respective  molecular  weights  nu  ,  are  used 
to  complete  the  density  calculation  by  a  simple  algebraic  formula. 

Implementation  of  the  ne  .•  model  produced  a  new  hybrid  SUBROUTINE 
DENSTY  to  replace  the  old  digital  subroutine  of  the  same  name.  For 
each  ^t  of  the  digital  solution  of  the  satellite  motion  equations, 
SUBROUTINE  DENSTY  is  called  to  compute  the  atmospheric  densi  v. 

In  the  subroutine,  T^  is  first  computed  on  the  digital  side  since 
the  needed  formulas  are  purely  algebraic.  T^  and  the  altitude,  r, 
are  then  transferred  to  the  analog  which  integrates  the  molecular 


density  equations.  Four  such  equations  are  solved  in  parallel  -- 
one  for  each  of  the  four  most  significant  gases  in  the  atmosphere: 

^2’  ^2*  anc*  mo^ecu^ar  densities  at  altitude  z  are 

transferred  back  to  the  digital  which  completes  the  density  calcu¬ 
lations. 

In  addition  a  density  check  subroutine  was  added  to  the  program  which 
when  called,  performs  a  dynamic  check  of  the  entire  hybrid  density 
model  to  confirm  that  it  is  operating  both  correctly  and  accurately. 

The  new  hybrid  orbiting  satellite  simulation  program  is  called 
SIMSAT  II. 

1.2.0  CALCULATION  OF  EXOSPHERIC  TEMPERATURE 

In  order  to  solve  for  the  number  densities  leading  to  the  air  density 
it  is  necessary  to  compute  the  exopsheric  temperature,  Tro  ;  i.e. 

Jie  temperature  at  the  outeimo  t  portion  of  the  atmosphere.  It  is  Tm 
which  is  modeled  as  a  function  of  all  the  parameters,  except  altitude 
which  influence  density.  The  effect  of  each  parameter  is  treated 
separately. 

1.2.1  VARIATION  WITH  SOLAR  ACTIVITv 

The  parameter  used  to  describe  variation  of  T  due  to  solar 

00 

activity  is  the  10.7-centimeter  solar  flux  which  is  monitored  by  the 
National  Research  Council  in  Ottawa,  Ontario.  Effects  of  the  slow 
11-year  cycle  variation  and  the  27-day  solar  rotation  are  both 
considered. 

1. 2. 1.1  VARIATION  WITH  THE  SOLAR  CYCLF 

Let  f  7  be  the  10,7-centimeter  solar  flux  in  units  of 
-2?.  ■'  2 

10  watts /m  /cycle/sec  and  T  the  night  time  global  minimum  value 

o 


of  ,  each  averaged  over  three  solar  rotations.  The  formula 
relating  these  two  quantities  for  quiet  geomagnetic  conditions  is: 

T  o  362  +  3.60  F,_ 
o  10.7 

I. 2. 1.2  VARIATION  WITHIN  ONE  SOLAR  ROTATION 

Let  F  ^  be  che  daily  mean  of  the  10.7-centimeter  solar  flux.  We 
can  correct  Tq  for  the  day-to-day  temperature  variation  due  to  the 
variation  within  one  solar  rotation.  The  formula: 

To'  "  fo  +  1'8  <*10.7  •  *10.7> 

I 

yields  Tq  ,  the  corrected  night  time  global  minimum  of  T^-  Values 

of  _  and  F,_  n  at  orbit  insertion  are  kept  constant  for  the 

id.  /  10. / 

lifetime  of  the  satellite. 


1.2.2  SEMI-ANNUAL  VARIATION 

The  semiannual  variation  is  described  by  the  function 
f (d)  =  ^0.37  +  0.14  sin  2rr  A|.*  ^  sin  4n 

where  d  is  the  number  of  days  elapsed  since  January  1  of  each  year. 


This  variation  is  used  to  produce  T  ,  the  night  time  global  minimum 
exospheric  temperature  according  to  the  formula 


T 

o 


The  combination  of  sine  terms  produce  two  unequal  pairs  of  maxima 
and  minima.  The  primary  and  secondary  maxima  occur  on  October  14  and 
April  20  respectively;  and  the  primary  and  secondary  minima  occur  on 
July  18  and  January  8  respectively. 
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Note  that  the  solar  and  semiannual  variations  can  be  combined  by 
substitution  into  the  single  formula: 

To  =  362  +  '•*  F10.7  +  (1'8  +  f(d))  F10.7 
1.2.3  DIURNAL  VARIATION 

The  distribution  of  exospheric  temperature  on  the  earth  is  such 
that  the  maximum  occurs  at  1400  hours  local  solar  time.  The 
atmosphere  bulges  out  in  the  bright  hemisphere  (the  dumal  bulge), 
producing  higher  densities  above  200  km.  The  center  of  the  bulge 
is  at  latitude  0  . 

D 

The  maximum  global  exospheric  temperature,  T  ,  and  minimum  global 
exospheric  temperature,  Tq,  are  related  by  the  formula 

T 

^  =  1  +  R 

o 

where  R  has  been  found  to  be  0.28. 

For  latitude  <f>  two  new  angles  as  defined: 

r  «  ;  |(MrJ  and  0  =  \ 

is  taken  to  be  zero;  i.e.  the  diurnal  bulge  is  modeled  as 
residing  on  the  equation.  Thus  n  =  Q  =  \  ( 0  |  . 

To  represent  the  effect  of  the  sun's  position,  an  angle  t  must  be 
computed : 

•A*  Vr 

t  <=  H  +  3  +  p  sin  (11  +  v  ) 

The  constants  used  in  the  above  formula  are  assigned  values  as  follows: 
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H  is  the  hour  angle  of  the  sun  measured  from  its  most  overhead 
position,  i.e.  1200  hours  local  solar  time.  The  relationship 
between  hour  angle  and  local  solar  time  is  illustrated  by  Figure  (1-1) 


Noon 
1200  LST 


0000  LST 
Midnight 


FIGURE  (1-1) 


The  computation  of  exospheric  temperature  is  completed  by  the  addition 

T  ■  T  +  M 

ro 

1.3.0  CALCULATION  OF  AIR  DENSITY 

Having  computed  the  exospheric  temperature,  the  deusity  for  any 
altitude  can  be  found  by  integrating  the  gas  equations.  The  solutions 
to  these  equations  produce  air  densities  that  are  considered  valid 


for  the  spring  and  fall  season.  For  the  summer  and  winter  seasons, 
modifications  must  be  made  to  complete  the  density  model. 

1.3.1  THE  DIFFUSION  EQUATIONS 


Diffusive  equilibrium  is  assumed  above  120  K  .  Using  the  ideal  gas 

m 

law  with  a  diffusion  term,  the  diffusion  equation  is  for  the  number 
density,  n^,  is  written  as: 


!I£  dz  .  dT 

K  T  T  (  1  +  ' 


(1) 


Here,  a  is  the  thermal-diffusion  factor,  K  is  the  Boltzmann  constant, 
and  g  is  the  acceleration  of  gravity,  nu  is  the  mo1eculat  weight  of 
the  constituent. 


T 

2a 

355 .0°K 

120 

,  „  11 

-3 

n  (No) 

SS 

4.0  x  10 

cm 

z 

n(02) 

*2 

7.5  x  1010 

-3 

cm 

n(0) 

- 

7.6  x  1010 

-3 

cm 

n(He) 

SC 

3.4  x  10^  cm  ^ 

s  is  a  function  of  T  and  is  found  by  the  formula: 


2 

s  =  0,0291  exp  (-q  /2)  in  units  1/km 


T  -  800 

on 

where  q  = - r - « 

750  +  1.722  X  10  (T  -  800) 

CO 


Because  great  variations  in  altitude  are  involved  with  these  calcu¬ 
lations,  the  acceleration  of  gravity,  g,  is  not  taken  as  being 
constant  but  is  found  from  the  inverse  square  formula: 


.8  = 


8  R 

o  o 


(R  +  z)‘ 


Where  g  is  the  value  of  g  at  the  earth's  surface  and  R  is  the 
°o  °  o 

earth's  radius,  taken  to  be  6371  km. 


Integrating  the  set  of  four  equations  (3)  from  120  km  to  altitude 

z  -/ields  the  needed  n.  . 

x 


The  density,  o,  is  then  found  by  the  simple  formula 


where  A  is  Avogadro's  number,  6.023  x  10  molecules /mole. 


1.3.2  THE  SUMMER  AND  WINTER  MODIFICATIONS 

It  is  found  that  the  densities  computed  from  the  diffusion  equations 
are  accurate  representations  of  the  behavior  during  spring  and  fall. 
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However,  modifications  must  be  made  for  the  other  two  seasons  because 
the  air  at  the  altitudes  considered  is  less  dense  in  summer  and  more 
dense  in  winter. 

The  deviations  from  the  spring-fall  models  are  illustrated  for  three 
values  of  by  the  chart  in  Figure  (1-2)  taken  from  Page  40  of  U.  S. 


PERCENT  DENSITY  DEPARTURES  FROM  SPRING/FALL 
MODEL8 


I  I CUR H  (I -2)-Densitv  departures  from  the  spring/fall 
models  for  summer  and  winter  with  three 
exospheric  temperatures. 
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Note  that  at  higher  altitudes,  the  winter  and  summer  densities 
merge  into  the  spring-fall  density.  The  lowest  altitude  that  the 
models  first  match,  Z  ,  is  a  function  of  T  and  are  taken  from 
Table  1  in  the  reference  and  are  listed  in  Table  (1-1). 


TABLE  (1-1) 


T 

CO 

CK 

2^  (Summer) 

km 

Z (Winter) 

km 

600 

195 

220 

700 

200 

225 

800 

210 

230 

900 

220 

235 

1000 

230 

240 

1100 

235 

240 

1300 

245 

245 

1500 

250 

250 

1700 

255 

255 

1900 

255 

255 

2100 

255 

255 

Below  Z^,  the  winter  density,  ow>  and  the  summer  density,  ,  are 
found  from  the  spring- fa  11  density,  n.  by  multiplying  by  the  factors 
found  in  the  formulas  below: 


where 


and 


W 

r 

S 

r 

-Ir 

D 


(1.4848  -  .4848,!,) 

(.7919  +  .2081, j,) 

tanh  1.75D  +  .0590 2 
z-120 

Z  -120 

m 


1.4.0  ORGANIZATION  OF  THE  HYBRID  PROGRAM 


The  air  density  model  described  was  incorporated  into  the  orbiting 
satellite  program  developed  previously  by  EAI  for  AFCRL  under 
Contract  Number  F1962S-67-C-0359.  The  all-digital  SUBROUTINE  DENSTY 
of  that  program,  which  modeled  density  only  as  a  function  of  altitude, 
was  replaced  by  a  new  hybrid  SUBROUTINE  DENSTY  to  implement  the 
more  realistic  model  described  in  this  report. 

The  day  of  the  year  and  Greenwich  Mean  Time  of  orbit  insertion  are 

read  in  on  cards,  as  are  both  the  daily  and  monthly  means  of  the 

10.7-centimeter  solar  flux.  Since  SUBROUTINE  INITAL  read  all  cards 

in  the  original  program,  these  new  read  statements  were  added  to 

that  subroutine.  SUBROUTINE  INITAL  also  calcuates  ^T  from  the 

value  of  K  read  in  on  a  ca^d.  Initialization  of  all  variables 
P 

used  in  the  density  calculation  is  done  by  this  subroutine. 

1.4.1  SUBROUTINE  DENSTY 

For  each  time  step  of  digital  integration  SUBROUTINE  DENSTY  is  called 
by  the  main  program.  The  subroutine  first  checks  if  satellite 
altitude  z  is  less  than  BRNOUT,  the  minimum  altitude  before  satellite 
burn-out.  If  it  is,  a  variable  IBNOUT  is  set  to  1,  and  the  main 
program  halts  the  simulation.  If  the  satellite  has  not  burned  out, 
the  subroutine  proceeds  with  the  density  calculation. 

First  the  day  of  the  year  is  calculated  from  the  initial  day  and 

the  total  elapsed  time.  The  exopheric  temperature  T  is  then  ealeu- 

00 

lated  using  all  the  algebraic  formulas  described  in  Section  2. 

This  is  used  to  compute  the  value  of  s.  Then  T^,  s,  and  the  satellite 
altitude  z  are  sent  to  the  analog  computer  through  digital-to-analog 
converters . 


A  control  line  is  reset  which  puts  the  analog  circuits  in  the 
operate  mode.  The  diffusion  equations  are  then  integrated  from 
120  km  to  z  .  (The  prime  is  used  to  distinguish  satellite 
altitude  from  the  independent  variable,  z,  on  the  analog  computer.) 
The  values  of  n^  for  Nj,  Oj,  0,  and  He  are  transferred  back 
to  the  digital  computer  by  analog-to-digital  converters  and  the 
analog  circuits  are  reset. 

The  spring-fall  density  is  then  calculated  as  a  weighted  sum  using 
the  molecular  weights  and  Avogadro's  number. 

The  proper  summer  or  winter  modification  is  multiplied  if  the 
season  is  not  spring  or  fall.  A  variable,  ICT,  indicates  the  number 
of  days  elaosed  in  the  present  season  of  the  year.  After  ICT  reaches 
91  the  subroutine  changes  the  season,  which  is  represented  by  the 
variable  ISN ;  the  seasons  summer,  fall,  winter  and  spring  being 
represented  by  the  values  1-4  respectively.  Note  that  a  seasonal 
change  (for  example  from  fall  to  winter)  would  add  an  "instantaneous" 
jump  in  density  if  the  modification  factor  were  applied  directly. 

In  order  to  smooth  out  this  jump  in  density  the  winter  and  summer 
modifications  are  weighted  by  a  smooth  curve  that  applies  them  such 
that  there  is  no  effect  at  the  beginning  of  the  season,  increases  to 
the  maximum  at  the  middle  of  the  season,  and  decreases  to  no  modifi¬ 
cation  again  at  the  end  of  the  season.  This  "modification  to  the 
modification"  is  described  by  a  second  order  algebraic  equation 
(parabola)  and  is  illustrated  in  Figure  (1-3). 

A  detailed  flow  chart  and  list  of  variable  definitions  for  SUBROUTINE 
DENSTY  are  given  in  Appendix  I A. 5. 


1007, 


Wr  or  Sr 


Day  of  Season 


Figure  (1-3)  Winter  and  Summer  Modification  Smoothing 
Function 


1.5.0  MODEL  VERIFICATION 


An  automatic  pot  setting  subroutine  and  static  check  subroutine 
have  been  incorporated  into  the  new  hybrid  orbiting  satellite 
program. 

In  addition,  a  complete  dynamic  check  of  the  density  model  is 
available  to  the  user. 


At  the  beginning  of  each  run,  the  operator  can  initiate  the  dynamic 


denrity  check  by  setting  a  console  switch.  The  analog  integrates 
the  diffusion  equations  up  to  a  number  of  different  altitudes  for 


each  of  several  values  of  T  .  The  log  n.  are  returned  to  the  digital 

CD  X 


and  the  density  computed  for  each  case. 


The  results  are  displayed  on  the  line  printer  in  the  same  format  as 
the  tables  at  the  end  of  the  U.  S.  Standard  Atmosphere  Supplements 
book.  Note  that  the  units  have  been  changed  for  the  printout  to 
agree  with  the  entries  in  the  tables.  This  allows  for  ready  com¬ 
parison  of  values  to  verify  that  the  model  is  operating  properly, 
and  if  not,  help  determine  the  source  of  the  problem. 


1.6.0  CONCLUSIONS  AND  RECOMMENDATIONS 


The  hybrid  air  density  model  implemented  is  a  great  improvement  over 
pure  altitude  functions  for  use  in  low  altitude  satellite  simulations. 
By  integrating  the  diffusion  equations  at  high  speed  on  the  analog 
computer,  the  digital  computer  can  still  integrate  the  satellite 
trajectory  equations  accurately  and  at  speeds  as  high  as  1000  time  real 
time. 


In  addition,  the  time  step  is  automatically  adjusted  to  be  smaller 
when  the  satellite  is  at  lower  altitudes.  Thus  maximum  accuracy  is 
achieved  at  the  portion  of  the  orbit  where  air  density  is  most  signi¬ 
ficant  . 


The  additional  parameters  included  insure  a  more  realistic  representa¬ 
tion  of  the  atmosphere  and  its  influence  on  satellite  motion.  However, 
to  further  improve  the  density  model  it  is  recommended  that  a  variable 

representation  be  included  for  the  parameters  F  ,  F.  _,  and  K  . 

10.7  iO.  /  P 

In  the  present  version  of  the  program,  the  values  of  these  three  parameters 
at  orbit  insertion  are  kept  constant  for  the  lifetime  of  the  satellite. 
Since  these  parameters  are  continuously  varying  an  improvement  would  be 
realized  by  replacing  them  with  time  dependent  functions  which  would 
better  represent  their  behavior. 
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IA1  AMPLITUDE  SCALING 

IA1.1  TEMPERATURE  EQUATION 

The  equation  for  temperature  variation  along  with  its  unsealed 
analog  circuit  diagram  are  shown  below: 


s 


The  scaling  table  is: 


Problem  Variable 
T 


Estimate  Maximum 
2500°K 

■1 


Computer  Variable 
T  T 


1/30  km 


The  scaled  equation  is  now  formed: 


dT 

dz 


-  s  (T  -  T  ) 

00 

-2500  ^fe]-^  H  \  [M  -  [4]} 

-  h  [boo]  ■  -°333  H  {  [2S00]  -  [4] } 

From  this  equation  the  scaled  diagram  is  drawn 


ADC 


The  scaling  table  is 


Problem  Variable  Estimated  Maximum 

z  1000  km 


The  scaled  equation  is  then  developed. 


Substituting  6371  for  Ro  yields 


The  scaled  diagram  is  thus: 


.6371 
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IA1.3 


DIFFUSION  EQUATIONS 


The  set  of  diffusion  equations  to  be  solved  are  repeated  here  for 
convenience : 


n.m.g  n.  (1+CA  s  (T-TI 
1  1  _i 

KT  +  T 


Because  of  scaling  considerations,  it  is  preferable  to  work  with 

In  n..  From  the  basic  relation  -r^—  In  n.  =  —  it  follows  that 
1  dn .  l  n . 

l  i 

dn^  =  n.  In  n.,  Substituting  in  the  above  diffusion  equations  yields, 
after  canceling  n. : 

l 


d  ( In  n  . ) 
dz 


,  (l  +  °0  s 

T 


(T  -  T  ) 

00 


It  is  this  final  set  of  equations  which  is  solved  on  the  analog 
computer. 

The  unsealed  diagram  that  describes  the  four  diffusion  equations 
is  shown  below: 
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Note  that  for  scaling  purposes  we  have  defined  the  outputs  of 
the  summer  and  divider  as  x  and  D  respectively 

In  order  to  be  consistant  with  the  table  entries  in  U.  S.  Standard 
Atmosphere  Supplements,  1966,  the  common  logarithm  is  substituted 
for  the  natural  logarithm  according  to  the  formula: 

In  n  =  2.3026  log  n 


The  unsealed  diffusion  equation  now  becomes: 

,  m  g  (1  +  a)  s  (T  -  T  ) 

2.3026  ±  (log  n)  =  -  ? - 


The  scaling  table  is: 


Problem  Variable 


Estimated  Maximum 


2500°K 


Computer  Variable 


1  u  -1 

30  km 


100  K/Km 


1  ,,  -1 
3  km 


log  n 


30 


The  scaled  equations  are  now  formed: 


n>ig 


-  ( 1  +  a)  s  (T  -  T  ) 


100 


x 

Too 


n^g 

Took 


miK  2500 
K  "  30 

,8333  (1  +a) 


[do-] 

[3°s] 

[™] 

N 


[*]■■(¥■¥“*->  W(fe] 

wife] -fell) 


r  T 

L2500. 


-2 -3026  (log  n)  =  D 
-30(2.3026,  fj  [l26-»]  .  >  [3D] 


t  fV]  '  -«0«26[»] 


From  these  equations  the  scaled  diagram  is  drawn: 


IA1.4  ALTITUDE  RAMP 

A  circuit  is  necessary  to  integrate  from  120  km  up  to  the  satellite 
altitude  z' . 

The  unsealed  equation  is  merely: 

dz  - 
dz  = 

But  z  has  a  scale  factor  of  1000  so  this  becomes 

dl  [1000]  *ool 

The  scaled  diagram  is: 


ALL  INT. 


Sense 


(This  pot  eliminated 
by  time  scaling) 
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IA2 .  "TIME"  SCALING  OF  THE  ANALOG  COMPUTER 

During  each  of  digital  integration,  the  digital  computer 
waits  for  the  analog  computer  to  integrate  the  diffusion 
equations  to  obtain  density.  Theiefore  to  keep  ^t  small  to 
insure  high  accuracy  it  is  advantageous  to  run  the  analog  computer 
as  fast  as  possible. 

A  value  for  p  of  10  ^  is  chosen  because  it  provides  the  best 
combination  of  speed  and  accuracy. 

According  to  the  formula 

t  =  pz 

This  means  that  z  will  integrate  1000  km  in  1  millisecond. 

1A3 .  ANALOG  COMPUTER  CIRCUIT  DIAGRAMS 

The  complete  scaled  diagrams  for  the  hybrid  density  model  are 
shown  in  figures  IA1  and  IA2.  Values  of  potentiometer  settings  and 
static  check  values  for  these  circuits  are  given  in  Tables  IA1  and 
IA2  respectively. 
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Figure  IA2  Analog  CU cults 


DAC 


m.g 

Calculation  of  — 

K 

m  =  m  grams /mole 
1  w 

where  m  is  molecular  weight  as  a  pure  number 
w  2 

g  =  .0098  km/sec 
o  20  2  2  o 

k  *  1.3805  x  10  gm  km  /sec:  K  molecule 


Using  above 


m.  g 

1  °  7  O' 

,rt23  molecules  °K 

K  ~  7-° 

7  o  x  iu  m  m  . 

w  mole  km 

We  can  simplify  the 
23 

10  molecules/mole 

units  by  dividing  by  Avogadro' 

which  yields 

s  number  6.02 

mi  8o 

-  1  1  a  m  -  --  - 

K 

X  i  lO  Ill  , 

w  km 

Using  this  formula, 

the  following  table  of  values 

is  formed 

Atmospheric 

M 

m.  g 

/i-mN 

Contituent 

w 

— ^  1  K./Km; 

N2 

28 

33.05 

0 

16 

18.88 

°2 

32 

37.76 

He 

4 

4.72 

30 


TABLE  (IA1 ) 

POTENTIOMETER  SETTINGS 
FOR  DENSITY  MODEL 


Pot  Number 

Parameter 

C600 

108  "120 

30 

C601 

.0004826 

P 

C602 

°2 

(scaling) 

C603 

mi  8o 

100K 

C610 

.8333  (1  +  a  ) 

C800 

108  "120 

30 

C801 

.0004826 

P 

C802 

He 

(scaling) 

C803 

m.  g 

1  0 

100K 

C810 

.8333  (1  +  a  ) 

TABLE  (IA1)  CONT'D 


I 


Setting 

.8625 

.4826 

.1200 

.3776 

.8333 

.7510 

.4826 
.  ?.200 

.0472 

.5166 


STATIC  CHECK 


Component 

Va  lue 

A002 

- . 4000 

A003 

.2580 

A012 

.8430 

R801 

-.2175 

P011 

-.0725 

POlO 

o 

CM 

• 

AOIO 

.1420 

AOll 

.2175 

P212 

.6371 

A212 

-.6491 

R804 

-.4213 

A213 

.4213 

P213 

-.4059 

M- 2F (231-R) 

.9634 

P411 

-.1200 

A411 

.1200 

P200 

-.8867 

A200 

.8867 

P201 

.2039 

M-3F  (231-R) 

.4225 

P202 

-.0600 

A203 

-.4996 

P203 

.3184 

P210 

.1812 

P400 

-.8627 

A400 

.8627 

P401 

.1482 

M-3B(231-R) 

.3070 

P402 

-.0436 

A403 

-.3631 

TABLE  (IA2) 
STATIC  CHECK  VALUES 
FOR  DENSITY  MODEL 


Component 

Va  lue 

P403 

.1819 

P410 

.1812 

P600 

-.8625 

A600 

.8625 

P601 

.2223 

M-3D  (231-R) 

.4606 

P6"2 

-.0654 

A6o: 

-.5450 

P603 

.3638 

P610 

.1812 

P800 

-.7510 

A800 

.7510 

P801 

.0642 

M-8F  (231-R) 

.1331 

P802 

-.0189 

A803 

-.1579 

P803 

.0455 

P810 

.1124 

TABLE  (IA2 )  CONT'D 
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IA4 . 


VARIABLE  DEFINITIONS  FOR  SUBROUTINE  DEN STY 


R 

R«J 

BRNOUT 

IBNOUT 

D 

IPD 

D0 

T 

ICT 

F 

T0 

GMT 

GMT0 

LAMN 

LSI 

IN 

ETA 

THETA 

TAU 

H 

A 

RATIO 

TNF 

TNFDLT 

TM 

Q 

s 

DA  (3) 
LOGN (4) 
RHO 


satellite  altitude 
orbital  radius 
radius  of  earth 

minimum  satellite  altitude  before  burn-out 

set  to  1  if  satellite  has  burned  out 

present  day  of  the  year 

integar  part  of  D 

day  of  year  of  orbit  insertion 

elapsed  time  in  seconds  from  orbit  insertion 

number  of  days  elapsed  in  present  season 

semiannual  variation  function 

night  time  global  minimum  exospheric  temperature 

present  Greenwich  Mean  Time 

Greenwich  Mean  Time  at  orbit  insertion 

longitude  of  satellite 

local  solar  time 

latitude  of  satellite 


V 

-  0  diurnal  bulge  angles 

TT 

* 

hour  angle  of  the  sun  (H  ) 

-  j  intermediate  values  in  diurnal 

-  )  variation  computation 

exospheric  temperature  (T  ) 

00 

geomagnetic  variation  (^T) 

(T  -  800) 

oo 


temperature  constants 


digital  to  analog  conversion  array 
analog  to  digital  conversion  array 
air  density  (o) 


ISN  -  season  of  the  year  indicator 

RNTNF  -  winter-summer  matching  altitude  pointer 

NTNF  -  integer  part  of  RNTNF 

ZMS(NTNF)  -  summer  altitude  matching  table 

ZMW(NTNF)  -  winter  altitude  matching  table 

ZM  -  matching  altitude, 

E  -  intermediate  value  in  summer-winter  modification 

PSI  -  summer-winter  modification  parameter 

SR  -  summer  modification  (Sr) 

WR  -  winter  modification  (Wr) 

A,B,C  -  constants  for  parabolic  smoothing  function 

SRMOD  -  modification  to  summer  modification 

WRMOD  -  modification  to  winter  modification 


IA6 


DIGITAL  LISTING  OF  SUBROUTINE  DEN STY 


SUBROUTINE  DENS TY ( R .LANIN »LN »RHO , I BNOUT , D ,DO ,GMTO * SFM ,SFD ♦ TNFDLT » 

1  I CT  »  I  SN  ♦  T  ) 

REAL  M,N»LN»i_AMN»LSTA 
REAL  LOGN I 4 ) 

INTEGER  G^TO,ZMS(l 1 ) »ZMW( 11 ) ,D,D0 
DIMENSION  DAO) 

COMMON  /CONSTS/  GO.RO.BRNOIJT 
DATA  PI/3.141593/ 

DATA  AR.M.N  /0. 28. 1 . 5 ,2 .5/ 

DATA  (ZMS( J) , J= 1,1 1 >/195* 200 .21 0,220 ,230,235 ,245.250*255*255,255/ 
DATA  (ZMW( J) ,J  = 1,1  I  1/2  20, 225, 23 0,235. 24 0,240*24 5*250*255, 255, 255/ 
DATA  BETA, P, GAMMA  /-0 . 78539 8 , 0. 20944 . 0 . 785 39 8/ 

E  30  =  1.0E30 

C  Z  IS  ALTITUDE  IN  KILOMETERS 
Z  =  .OOi#(R-RO) 

IF  ( Z.LT.BRNOUT  )  GO  TO  30 
I  PD  =  D 

ID  -  DO  ^  T/86400.0 
0  =  MODI  ID, 36 5) 

IF  (D-IPD.EQ.l)  ICT  =  ICT  +  1 

F  =  I  0.37-t-O.  1 4*S I N  I  2.0*PI  *(  D-15I  )  /  365  )  )  *S  I  N  (  4«0*PI#(D-59)  /  365  ) 

TNFO  =  362.0  +  1.8*SFD  +  (1.8+F)*SFM 

GMT  =  GMT  0  +  T/ 36. 0 

LSTA  =  GMT  +  LAMN*2400.0/(2.0*PI ) 

LST  =  AMODILSTA  ,2400.0 )  0.5 

ETA  =  0 .5  *ABS ( L  N ) 

THETA  =  ETA 
C  H  IS  H* 

H  =  ( (LST-1200) /1200)*PI 

TAU  =  H  +  BETA  +  P*5 I N ( H+GAMMA ) 

IF  (TAU.GT.PI)  TAU  =  TAU  -2.0*PI 
IF  (TAU.LT.-PI)  TAU  =  TAU  +  2.0*PI 

A  =  AR* ( (  I  COS (ETA)  ) **M  -  ( S I N I  THE T A ) ) **M ) / I 1  +  AR* ( S I N ( THETA ) ) **M ) ) 
RATIO  =  (  H-A*(C05( TAU/2.0) )**N)*( 1+AR*(SIN(THETA) )**M) 

INF  =  T NF 0*RAT I  0 

TMF  =  TNF  +  TNFDLT 

IF  ( TNF. LT. 600.0)  TMF  =  600.0 

IF  ( TNF. GT. 2100.0)  TNF  =  2100.0 

IF  ( TNF. GT. 1100.0)  GO  TO  5 

RNTNF  =  (TNF-500. 0/100.0 

GO  TO  6 

5  RNTNF  =  ( TNF+1G0.0 ) /200.0 

6  CONTINUE 

TM  =  TNF  -  800.0 
q  =  TM/I750.0  +  .0001722*TM*TM) 

S  =  ,8730»EXPI -Q^0/2.0) 

D  A ; i )  =  -s 

DA<2)  =  TNF/ 2500.0 

DAO)  =  Z/1000.0 

CALL  QOAJMOIOO  ,D  A  ,  I  ERR  ,  JCH  AN  ) 


ASSEMBL 

TSL 

=  •6041  ♦ ,  2 

RESET  SENSE  LINE  1 

SFL 

=  *6041  ,  ,2 

RESET  CL  1  OPERATE 

STIJ  TSL 

=  *6041  ,  ,2 

TEST  SL1 

'  *1  IfPWSWlllW 
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STU 


JNZ 

FORTRAN 

CALL  QADCV0(0*4*LOGN*IERR*JCHAN) 

/  SFL  =’6061  *  *2  SET  CLl  IC 

RHO  =  4.6'315t-35*(  F30)**LOGN(  1  )  '+  2 .656  1 E-35*  (  E30  )  **LOGN  (  2  ) 
1  +  5.3122t-35*(E30)**LOGN(3)  +  0.66453E-35* ( E30 ) **LOGN ( 4 ) 

C  SEASON  SUMMER  FALL  WINTER  SPRING 

C  ISN  1  2  3  4 

C  ICT  IS  HOW  MAN',  DAYS  INTO  THE  SEASON 
IF  ( ICT.LT.92)  GO  TO  9 
ISN  =  MQD( ISM *4)  +  1 
ICT  =  1 

9  GO  TO  (70*71*72*73) *ISN 

70  MTNF  =  RNTNF 
ZM  =  ZMS(NTNF) 

IF  (Z.GE.ZM)  GO  TO  73 
E  =  ( Z-l 20 • 0 ) / ( ZM- 120.0) 

PSI  =  T  ANH ( 1 »75*E)  +  0.059*E*E 
SR  =  0.7919  +  0 . 20  8  1  *P  S I 
A  =  .000493827*1 1,0-SR ) 

B  =  -92. 0*A 
C  =  1.0-A-B 

SRMOD  =  A  *  I C  T  *  I C  T  +  B*ICT  +  C 

RHO  =  SRMOD*RHO 

RETURN 

71  CONTINUE 
73  RETURN 

72  NTNF  =  RNTNF 
ZM  =  ZMW(NTNF) 

r  F  (Z.GE.ZM)  GO  TO  73 
E  =  (Z-120.0)/(ZM-120.0> 

PSI  =  T  ANH ( 1 • 75  *E )  +  0.O59*E*E 
WR  =  1.4848  -  0.4848*PSI 
A  =  .000493827*( 1.0-WR ) 

B  =  -92 • 0*A 
C  =  1.0-A-R 

WRMOO  =  A# i CT* I CT  +  B* I CT  +  C 

RHO  =  WRMOD*RHO 

RETURN 

30  IBNOUT  *  1 
RETURN 
END 
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OPERATING  INSTRUCTIONS  FOR  SIMSAT  II 


The  orbiting  satellite  simulation  program  with  the  hybrid  air 
density  model  has  been  given  the  name  Simsat  II. 


For  the  sake  of  completeness,  the  operating  instructions,  although 
an  independent  self-contained  document,  are  included  in  this  section. 


OPERATING  INSTRUCTIONS 


FOR 

S  I  M  S  A  T  II 

(Satellite  Trajectory  Simulation  Program) 


Prepared  by 

Stuart  B.  Mindlin 

E1ECTRONIC  ASSOCIATES,  INC. 
Princeton  Computation  Center 
Princeton,  New  Jersey  08540 


Prepared  for 

AIR  FORCE  CAMBRIDGE  RESEARCH  LABS. 

United  States  Air  Force 
Bedford,  Massachusetts  01730 


October  8,  1969 


SIMSAT  II  is  a  satellite  trajectory  simulation  program  written  for  the  8900  hybrid 
computer.  This  program  was  designed  for  use  in  evaluating  the  effects  of  various 
perturbative  forces,  such  as  atmospheric  drag,  solar  radiation  pressure,  and 
gravitational  anomalies,  upon  the  trajectory  and  lifetime  of  an  orbiting  satellite. 
SIMSAT  II  is  identical  to  SIMSAT  I  with  the  addition  of  a  more  sophisticated 
density  model  as  discussed  in  U.  S,  Standard  Atmosphere  Supplements,  1966.  Inte¬ 
gration  of  equations  to  calculate  atmospheric  density  and,  in  turn,  drag  force 
is  done  on  the  8800  anaiog  portion  of  the  hybrid  system.  In  addition  variables 
of  interest  are  transferred  to  the  8800  analog  computer  where  they  may  be  ob¬ 
served  and  plotted.  To  operate  SIMSAT  II,  the  following  steps  must  be  followed. 

1)  Patch  a  pair  of  8800  analog  and  logic  panels  as  shown  in  Figures  TAi,  IA2, 
and  IA3.  These  are  the  equations  for  molecular  density  of  the  major 
atmospheric  constituents.  The  additional  patching  allows  the  user  to  plot 
latitude  versus  longitude  on  the  x-y  plotter  and  also  to  control  the  transfer 
of  program  variables  to  the  analog  computer. 

2)  Punch  a  set  of  data  cards  which  specify  initial  position  and  velocity  of 
the  satellite,  and  the  scale  factors  to  be  used  in  digital/analog  transfers, 
as  well  as  initial  data  needed  for  the  density  calculation.  (See  Table  2, 
section  7.) 

3)  Mount  the  8800  patch  panels.  Switch  1013  off. 

4)  Load  and  execute  the  digital  program. 

5)  Control  of  the  running  program  may  be  maintained  from  either  the  analog  or 
digital  consoles  as  follows: 

a)  If  console  register  switch  8  is  on,  control  is  from  the  digital  console, 
otherwise  it  is  from  the  analog  console, 

b)  If  control  is  from  the  digital  console  (C8  on),  then  the  mode  of  the 
simulation  is  determined  by  the  following  table. 


Mode 

Console  Switch 

1 

Console  Switch 

2 

IC 

ON 

OFF 

IC 

ON 

ON 

HOLD 

OFF 

OFF 

OP 

OFF 

ON 

If  it  is  desired  to  print  out  the  values  of  the  program  variables  on 
the  6400  line  printer,  turn  on  console  register  switch  7.  All  program 
variables  will  have  their  values  printed  out  once,  and  the  console 
switch  will  be  reset.  This  may  be  done  whether  the  simulation  is  under 


c) 


analog  or  digital  control,  but  only  when  in  OP  or  HOLD  modes, 

d)  If  it  is  desired  to  read  in  new  initial  conditions  for  the  satellite, 
this  may  be  done  by  pressing  console  register  switch  10,  when  the 
simulation  is  in  IC  mode,  under  either  analog  or  digital  control, 

e)  If  console  register  switch  9  is  set,  the  program  will  go  into  HOLD 
mode  after  each  complete  orbit  of  the  satellite.  This  is  useful  in 
checking  out  the  change  in  orbital  parameters  from  one  orbit  to  the 
next . 

f)  Which  program  variables  are  transferred  to  the  analog  computer  is 
controlled  by  four  logic  function  switches  on  the  8800  console,  as 


follows : 

DAC  CHANNEL 

FUNCTION  SWITCH  Oil 

LEFT  OR  CENTER 

FUNCTION  SWITCH  Oil  RIGHT 

3 

radius 

radial  perturbation 

4 

longitude  in  orbital  plane 

perturbation  in  orbital  longitude 

5 

velocity  normal  to  radius 

perturbation  in  tangential  velocity 

6 

velocity  in  radial  direction 

perturbation  in  radial  velocity 

7 

geographical  longitude 

Keplerian  radius 

8 

geographical  latitude 

Keplerian  orbital  longitude 

9 

altitude 

Keplerian  tangential  velocity 

10 

air  density 

Keplerian  radial  velocity 

11 

x  component  of  drag  force 

x  component  of  drag  force 

12 

y  component  of  drag  force 

y  component  of  drag  force 

13 

x  component  of  velocity  See 

tangential  velocity 

14 

y  component  of  velocity  Table 

zero 

15 

z  component  of  velocity  1 

radial  velocity 

As  outlined  in  the  report  AFCRL  69-0123 
used  in  the  simulation.  The  satellite': 

several  coordinate  frames  are 
s  velocity  is  resolved  into  each 

of  these  coordinate  systems,  and  the  x,y,  and  z  components  of  any  of 
these  velocity  representations  may  be  transferred  on  DAC  channels  13,  14, 
and  13,  respectively.  Which  velocity  components  are  transferred  is 
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controlled  by  logic  switches  1012,  1013.  and  1211  a3  follows: 


TABLE  1 


Switch  1012 

j  Switch  1013 

Switch  1211 

Velocity  Transferred 

OFF 

OFF 

OFF 

H  frame  components 

OFF 

OFF 

ON 

H  frame  components 

OFF 

'  ON 

OFF 

.  Modified  Euler  frame  components  i 

OFF 

ON 

ON 

1  I 

J  Orbital  frame  components  j 

t 

ON 

OFF 

OFF 

Inertial  frame  components 

ON 

OFF 

ON 

Navigational  frame  components 

ON 

ON 

OFF 

Elite**  frame  components 

.  -  .  -  ---  -  .j 

ON 

>  on 

ON 

Euler  frame  components  J 

7)  The  format  of  the  data  cards  required  by  SIMSAT  II  is  described  below. 

The  first  four  cards  contain  the  scale  factors  to  be  used  in  the 
D/A  transfer  of  program  variables.  Cards  1  and  2  contain  the  scale 
factors  to  be  used  when  logic  switch  1012  is  on,  and  cards  3  and  4 
contain  the  scale  factors  to  be  used  when  switch  1012  is  off.  Each 
scale  factor  is  punched  in  10  column  floating  point  format.  The 
fifth  card  contains  gQ,  the  gravitational  constant  in  columns  1-10, 
R0,  the  radius  of  the  earth,  in  columns  11-20,  and  %in»  the  altitude 
below  which  the  satellite  is  considered  to  have  burnt  out,  in 
columns  21-30.  During  execution  of  SIMSAT  II,  if  the  altitude  of  the 
satellite  goes  below  Hmin»  the  run  Is  terminated  and  the  time  at 
which  burnout  occurred  is  printed  out.  The  last  three  values  on 
the  fifth  card  are  for  use  in  the  density  calculation.  The  3-hour 
geomagnetic  planetary  index,  Kp,  is  found  in  columns  31-40.  Coluims 
41-30  and  51-60  contain  respectively  the  monthly  and  daily  means  of 
the  10.7  cm.  solar  flux.  (See  U.  S.  Standard  Atmosphere  Supplements, 
1966,  pg.  47.)  Card  number  6  contains  either  "YES'*  or  "NO"  in 
columns  J-3,  to  indicate  if  the  satellite's  initial  position  and 
velocity,  orbital  insertion  time  and  day  are  to  be  read  in  from 
cards,  or  typed  in  on  the  console  typewriter.  NO  means  the  data 
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is  to  be  read  in  from  cards,  in  which  case  three  more  data  cards  are 
required  which  contain  the  following  initial  data,  the  first  two  of 
which  are  in  10  column  floating  point  format:  Badius  in  meters, 
longitude  in  degrees,  latitude  in  degrees,  velocity  in  meters/sec., 
azimuth  in  degrees,  flight  path  angle  in  degrees,  mass  in  kg.,  area 
in  square  meters,  drag  coefficient,  and  time  scale  factor  (factor 
by  which  simulation  is  to  be  sped  up  over  real  time,  normally 
(1000.0).  The  third  card  contains  the  day  of  the  year  in  columns 
5-7  and  Greenwich  Mean  Time  at  orbit  insertion  in  columns  11-14. 

The  layout  of  all  the  data  cards  is  summarized  in  Table  2. 


TABLE  2 

Scale  Factors  Should  be  the  Reciprocals 


Card 

Columns 

" 

of  the  1 

Maximum  Valve 

1 

1-10 

Scale 

factor 

for 

DAC3, 

switch 

Oil 

right. 

1 

11-20 

Scale 

factor 

for 

DAC4, 

switch 

Oil 

right. 

1 

21-30 

Scale 

factor 

for 

DAC5, 

switch 

Oil 

right. 

1 

31-40 

Sea  le 

factor 

for  DAC6, 

switch 

Oil 

right . 

1 

41-50 

Scale 

factor 

for 

DAC7, 

switch 

Oil 

right . 

1 

51-60 

Sc3  ie 

factor 

for  DAC8, 

switch 

Oil 

right. 

1 

61-70 

Scale 

factor 

for 

DAC9, 

switch 

Oil 

right. 

1 

71-80 

Sea  le 

factor 

for 

DAC10 

,  switch 

Oil 

right . 

2 

1-10 

Scale 

factor 

for 

DAC11 

,  switch 

Oil 

right. 

2 

11-20 

Sea  le 

factor 

for  DAC12 

,  switch 

Oil 

right. 

2 

21-30 

Sea  le 

factor 

for 

DAC13 

,  switch 

Oil 

right. 

2 

31-40 

Sea  le 

factor 

for 

DAC14 

,  switch 

Oil 

right. 

2 

41-50 

Sea  le 

factor 

for 

DAC15 

,  switch 

Oil 

right. 

2 

51-60 

(Blank) 

2 

61-70 

(Blank) 

2 

71-30 

^Blankl 

3 

1-10 

Sea  le 

factor 

for 

DAC3, 

switch 

Oil 

left  or 

center . 

3 

11-20 

Sea  le 

factor 

for 

DAC4, 

switch 

Oil 

left  or 

center. 

3 

21-30 

Sea  le 

factor 

for 

DAC5 , 

switch 

4 

left  or 

center. 

3 

31-40 

Sea  le 

factor 

for 

DAC6, 

switch 

on 

left  or 

center. 

3 

41-50 

Sea  le 

factor 

for 

DAC  7 , 

switch 

Oil 

left  or 

center. 

3 

51-60 

Sea  le 

factor 

for 

DAC  8 , 

switch 

ou 

left  or 

center. 
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TABLE  2  (continued) 


Card 

Columns 

Data 

3 

61-70 

Scale  factor  for  DAC9,  switch  Oil  left  or  center. 

3 

71-80 

Scale  factor  for  DAC10,  switch  Oil  left  or  center. 

4 

1-10 

Scale  factor  for  DAC11,  switch  Oil  left  or  center. 

4 

11-20 

Scale  factor  for  DAC12,  switch  Oil  left  or  center. 

4 

21-30 

Scale  factor  for  DAC13,  switch  Oil  left  or  center. 

4 

31-40 

Scale  factor  for  DAC14,  switch  Oil  left  or  center. 

4 

41-50 

Scale  factor  for  DAC15,  switch  Oil  left  or  center. 

4 

51-60 

(Blank) 

4 

61-70 

(Blank) 

4 

71-80 

(Blank) 

5 

1-10 

gQ,  gravitational  const,  (m/sec. ^). 

5 

U-20 

Rq,  radius  of  earth  (meters). 

5 

21-30 

burnout  altitude  (kilometers). 

5 

31-40 

Kp,  geomagnetic  planetary  index. 

5 

41-50 

Monthly  mean  of  10.7  cm  solar  flux. 

5 

51-60 

Daily  mean  of  10.7  cm  solar  flux. 

6 

1-3 

"YES"  or  "NO",  input  is/is  not  from  typewriter. 

7 

1-10 

Initial  radius  (meters). 

7 

11-20 

Initial  longitude  (degrees  east  of  Greenwich). 

7 

21-30 

Initial  latitude  (degrees  north  of  equator). 

7 

31-40 

Initial  velocity  (meters/sec.). 

7 

41-50 

Initial  azimuth  (degrees  clockwise  from  north). 

7 

51-60 

Initial  flight  path  angle  (degrees  away  from  earth) 

7 

61-70 

Satellite  mass  (kilograms). 

7 

71-80 

Satellite  cross-sectional  area  (square  meters). 

8 

1-10 

Drag  coefficient. 

8 

11-20 

Time  scale  factor. 

9 

5-7 

Day  of  the  year. 

9 

11-14 

Greenwich  Mean  Time. 

If  card  no.  6  is  YES,  indicating  typewriter  input,  the  user  should 
follow  the  directions  typed  out,  and  input  the  initial  data  as  it 
is  requested. 

8.  After  new  initial  conditions  are  read  the  system  will  type  that 
tie  operator  may  set  a  certain  console  switch  for  a  density  check. 

If  the  operator  sets  this  switch  the  system  will  automatically 
perform  a  dynamic  check  of  the  entire  hybrid  air  density  model. 

The  results  of  the  check  is  output  on  the  lineprinter  in  the  form 
of  the  tables  of  SPRING-FALL  density  in  the  U,  S.  Standard 
Atmosphere  Supplements,  1966 . 

Comparison  of  these  values  will  confirm  proper  operation  of  the 
model,  or  if  there  are  difficulties,  help  locate  the  source. 
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PROPAGATION  OF  ELECTROMAGNETIC 
ENERGY  IN  THE  IONOSPHERE 


II. 1.0 


INTRODUCTION 


This  section  of  the  report  describes  modifications  which  have  been 
made  to  the  basic  hybrid  ray  tracing  prograr..  that  was  developed 
under  a  previous  contract  (No.  F  19628-67 -C-0359) .  The  results 
and  derivations  of  the  ray  tracing  equations  were  included  in  the 
final  report  on  that  contract:  AFCRL  Report  69-0123,  "Hybrid  Com¬ 
puter  Applications  to  Mathematical  Models  of  Physical  Systems," 
February  1969. 

The  modifications  to  the  program  were  as  follows: 

(1)  Addition  of  ele  .tron  collisions. 

(2)  Addition  of  an  automatic  covergence  scheme 
based  on  the  range  between  transmitter  and 
receiver. 

(3)  Time  history  storage  of  pertinent  variables  for 
the  optimum  ray  path  when  determined  by  the 
covergence  scheme. 

1 1. 2.0  DESCRIPTION  OF  ELECTRON  COLLISIONS 

The  ray  tracing  equations  solved  to  determine  the  path  of  electro¬ 
magnetic  energy  transmitted  through  the  ionosphere  are  describ  '  in 
AFCRL  Report  69-0123,  Section  III.  In  these  equations  the  expression 
for  the  refractive  index  neglected  the  effect  of  electron  collisions. 
Under  the  current  contract  the  Appleton -Hartree  formula  for  the 
refractive  index  was  modified  to  include  the  general  case  of  electron 
collisions.  This  formula  in  genera1  form  is: 
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whe  * o , 

u  =  refractive  index 

,  2 
t.. 


fN  =  8.98  x  10  x  N‘,  the  plasma  frequency 

f  =  wave  frequency 
N  =  electron  density  (electrons /cc) 

Y  =  normalized  magnitude  of  the  earth's  magnetic  field  vector 

1  eB 

2-Tj;  m 

ip  ~  angle  between  "y,  earth  magnetic  field  vector,  and  the 
wave  normal 
Yt  -  Y  sin  « 

Yl  *  Y  cos  ,ir 

z  =  electron  collision  frequency 
+  p  an  ordinary  ray 
-  an  extra  ordinary  ray 

1 1. 3.0  DESCRIPTION  OF  CONVERGENCE  SCHEME 


The  convergence  scheme  used  is  basically  a  gradient  method  which 
corrects  the  firing  angle  («)  of  the  transmitter  based  on  the 
longitude  and  altitude  of  the  ray  path  upon  arrival  at  the  latitude 
of  \FCRL.  At  the  outset  it  was  thought  desirable  to  correct  not  only 
o,  But  also  the  azimuth  angle  (3).  However,  it  was  later  decided 
to  maintain  3  constant  since  magnetic  field  effects  were  not  thought 
to  be  of  significant  magnitude  as  to  cause  large  deviation  ir.  3 
along  the  path  of  the  ray. 


By  way  of  introduction  to  the  scheme  used,  it  should  first  be 
pointed  out  that  for  a  given  electron  density  distribution,  the 
approximate  range  of  firing  angle  is  not  in  general  known.  For 
this  reason,  it  was  thought  best  to  choose  three  firing  angles 
separated  sufficiently  far  apart  so  as  to  determine  the  approxi 
mate  area  in  which  to  advance  (see  flow  chart  and  convergence 
scheme  flow  for  specifics  of  program). 

The  first  three  paths  for  a  given  electron  density  profile  are 
shown  qualitatively  below.  In  this  instance,  the  path  which 
is  closest  to  the  desired  path  is  curve  number  3. 


FIGURE  (II-l) 
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CONVERGENCE  SCHEME  MATHEMATICS 


If  a*,  p*  are  defined  to  be  the  firing  and  azimuth  angle, 
i espectively ,  which  gave  the  minimum  altitude  at  AFCRL  latitude, 
we  write 

(1)  D(«*,  p*)  *  H(a*,  p*)  +  A  (9-9t) 

where  A  is  a  constant  to  hi.  defined  which  will  dictate  the  necessary 
tolerance  on  0  (longitude),  H  (a*f  P*)  is  the  final  height,  and 
0T  =  longitude  (target)  of  AFCRL. 

D  is  the  quantity  which  is  monitored  for  the  optimum  ray  path,  i.e. 
if  we  require  a  1  km  tolerance  on  the  final  altitude  H  (Cl,  P)  and 

a  .1  tolerance  on  (0-q  )  then  equation  (1)  becomes 

T 

(2)  D(a*,  p*)  =  l  +  l(.l)  =1.1 

so  that  whan  a  ray  path  is  found  whose  final  conditions  satis 'ey  (2), 
then  this  is  the  optimum  ray.  As  it  turns  out  magnetic  field 
deflections  are  slight  and  as  a  result  the  boundary  condition  of 
equation  2)  is  effectvely  a  function  of  final  altitude  only,  i.e. 
stopping  the  hybrid  program  at  the  latitude  of  AFCRL  also  means 
stopping  at  the  correct  longitude. 

Of  course,  in  general  one  of  the  first  three  rays  will  give  us 
only  an  approximation  of  the  general  area  in  which  to  advance. 
Therefore,  since  equation  (1),  (2)  will  most  often  not  be  satisfied, 
the  automatic  gradient  scheme  takes  over  to  choose  a  new  alpha  and 
beta  increment  as  follows:  To  start  the  pror _ ss  O  and  P  are 
incremented  (for  the  general  case;  only  a  was  incremented,  6  was  held 
constant)  by  a  constant  (1°)  to  determine-  to  slope.  First  increment 
aipha  and  monitor  the  value  of  D.  he  then  have 

(3)  D(«*+A«,P*)  .  H(a*+A«,p*)  +  A  (0(a*+A«,p*)  -9t) 


5 


and  as  an  approximation, 


(4) 


_  p(g*  +  £&t p*)  -  p(a*,p*j 
W  "  Aa 


Now  increment  P  and  write, 


(5)  D(a*,p*  +  Ap)  *  H(a*,p*  +  AP)  +  A(9<2*,0*  +  AP)  -  0T) 


3D  _  r)(ct*,P*  +  AP)  -  D(a*,P*) 


using  this  information  we  must  now  choose  the  next  <2,  P  pai. 


The  nugntidue  of  the  gradient  of  D  is  defined  as 

(7,  +f^ 


(8)  \»  m 


Therefore  we  may  write 


(9)  Aa 
K  ’  u  NEW 


H(«*  +  Aa.p*)  3D 


(10)  iP  -  ±jm  m. 

NEW  ,-m2  JP 


which  reduces  to  (when  =  0) , 


(U) 
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and 


<12>  *w-° 

Using  the  Aa  the  hybrid  program  starts  anew  and  the  same  process 
is  repeated  until  convergence  criteria  (equations  1,  2)  are  met. 

Once  the  optimum  ray  path  is  reached  it  is  repeated  so  that  the 
time  history  array  ra}  be  stored  and  printed  at  the  end.  This  is 
the  end  of  one  run  for  a  given  electron  density  profile. 

1 1. 3. 2  CONVERGENCE  FCLEME  FLOW 

As  mentioned  in  the  introductory  remarks  to  this  section,  the 
first  phase  of  the  convergence  scheme  flow  is  to  pick  three 
firing  angles  and  determine  the  ray  path  for  each.  These  angles 
as  may  be  seen  from  the  flowchart  Apendix  A3  are  12.5,  25,  37.5 
degrees.  Referring  to  the  system  flowchart,  at  statement  615 
we  begin  the  scheme.  After  first  having  chosen  Ot,  set  the 
necessary  pots,  and  initialized  the  program,  we  enter  the  hybrid 
loop  at  statement  number  70.  TCLl  is  the  Fortran  constant  which 
specifies  the  tolerance  on  the  latitude  as  a  stopping  condition, 
i.e.  once  DIFF,  which  is  the  difference  between  the  present 
latitude  and  the  final  latitude,  is  less  than  or  equal  to  T0L1 
we  fall  out  of  the  hybrid  loop  and  call  subroutine  QHOLD0.  This 
subroutine  (which  is  part  of  the  hybrid  run  time  library)  places 
the  analog  computer  in  the  HOLD  mode. 

The  next  decision  block  is  necessary  for  the  time  history  storag< 
END  is  a  logical  variable  which  is  set  true  upon  completion  of 
the  run  which  satifies  the  criterion  D(Of*,  P*)  <1.1  as  expressed 
by  equation  (2)  of  section  (3.1).  Once  END  becomes  true  we  repeat 
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the  run  so  that  we  may  store  the  pertinent  time  history 
variables.  The  criterion  expressed  above  is  stated  in  the  Fortran 
program  as  DIFFD  <;  T0L2;  hence  the  decision  block  for  DIFFD. 

If  DIFFD  is  greater  than  TOL  2  and  the  first  three  firing  angles 
have  been  completed,  we  incr^nent  0!  (statement  number  651)  and 
proceed  through  the  program  again.  The  remaind2r  of  the  flowchart 
follows  exactly  the  mathematical  description  in  the  last  section. 
This  procedure  begins  at  statement  numbei  1701.  The  following 
equivalences  are  helpful  in  comparing  the  FORTRAN  program  to  the 
matl'e-iatics : 

DADA  =  D(a*,  ft*) 

PDPA  =  9D 

3a 

PDPB  »  M 

SP 

2  2  2 
GRADD2  =  VD  =  M 

**  ap 

GRADD  *  vD 
DALPH  •= 

DBETA  =  AP 
AORIG  »  a  original 
BORIG  =  P  original 

Note  the  presence  of  he  logical  variable  NOTBE  on  the  third  page 
of  the  flowchart.  This  variable  acts  as  an  indicator  in  determining 
if  P  is  to  be  corrected  along  with  O'.  As  mentioned  previously,  3 
optimization  was  not  desired,  thus  NOTBE  was  false  during  all  runs. 

I I. 4.0  RESULTS 

If  the  ray  paths  are  not  convergent  the  automatic  scheme  can  do 
nothing.  In  this  case  there  is  no  time  history  stored.  In 
running  the  various  twelve  electron  density  profiles  eight  of 


twelve  profiles  would  not  allow  convergence  to  occur  at  the  given 
altitude,  latitude  and  longitude  of  AFCRL.  Essentially  what 
happens  is  that  the  program  chooses  the  next  angle  based  on  the 
previous  calculation  and  finds  that  instead  of  decreasing  the  final 
altitude  it  increases  it,  i.e.  the  final  altitude  as  a  function 
of  firing  angle  «  reaches,  at  least,  a  local  minimum.  This  event 
is  described  qualitatively  by  Figure  (II-2). 


In  this  case  the  program  tries  again  to  decrease  the  height  by 
going  in  the  other  direction  and  this  process  becomes  endless. 

On  most  of  the  rays  which  did  not  converge  the  final  altitude  was 
in  the  range  of  100-150  km. 


These  results  were  verified  by  the  conventional  ray  tracing  program 
with  no  convergence  scheme. 


Based  on  the  fact  that  the  automatic  scheme  showed  non-convergence 
for  eight  of  twelve  electron  density  profiles  it  was  desired  by 
AFCRL  to  have  the  running  procedure  for  the  conventional  hybrid 
ray  tracing  program.  The  procedures  supplied  for  both  programs 
are  included  in  Appendix  II  . 
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II.  5.0 


CONCLUSIONS  AND  RECOMMENDATIONS 


Based  on  the  results  of  the  convergence  scheme,  it  might  be 
desirable  to  first  make  small  modification  to  the  automatic 
program  which  would  all- w  the  program  to  pick  say  the  ray  path 
which  gave  the  absolute  minimum  altitude  at  the  latitude  and 
longitude  of  AFCRL,  for  a  given  profile.  Once  this  :‘.s  completed, 
it  would  be  informative  to  ±-  treduce  perturbations  into  the  electron 
density  distribution  during  a  transmission  for  the  ray  path 
determined  above.  The  results  of  this  might  be  used  for  comparison 
with  ionosonde  data.  For  example  what  perturbations  (or  random 
disturbance  pattern)  introduced  into  the  program  correlate  most 
nearly  with  ionosonde  data  regarding  intensity,  dispersion  etc. 

The  results  of  study  of  this  nature  might  be  very  useful  in  pre¬ 
dicting  at  least  qualitative  changes  in  electron  density  distributions. 


APPENDIX  II 


IIA1. 


ANALOG  SCHEMATICS 


Figures  (IIA1.)  and  (I1A2.)  are  the  analog  diagrams  of  the  automatic  system. 


(IIA1.)  Main  Analog  Program 


IIA2. 


COMPONENT  ASSIGNMENTS 


Pages  65  through  69  are  the  pot  and  amp.  assignment  sheets  plus 
static  test  voltiges.  Page  70  is  the  switch  assignment  sheet  for 
the  sense  lines  on  the  8800. 
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Pot  Number 

Parameter 

Setting 

C203 

Scaling 

.1111 

C210 

18 

4 

10  0 

.1800 

C310 

Ho  (IC) 
io3 

.0500 

C410 

Scaling 

.5000 

C601 

57.3 

2xl05  0 

.0286 

C810 

57.3 

4xl05  0 

.0143 

C812 

ina  1 

1000 

.6070 

C001 

1 

100  0 

1.0000 

CO  10 

1 

1000  0 

.1000 

C603 

2 

1000  0 

.2000 

C513 

Scaling 

.2000 

C612 

Sea  ling 

.2000 

Pot  Number 


Parameter 


Setting 


C213 

Sea  ling 

.0010 

C312 

H, 

*Scaling  for  ht  .  f 

1000 

.0100 

C311 

Scaling 

.  1000 

C701 

_2£_ 

.2500 

200 

C’910 

4>o 

.1000 

4 00~ 

*Hf  iS  the  hei§hc  at  which  there  is  a  scale  change  on  H  for  better 


resolution  H  =  10  km. 


...  . 


STATIC  CHECK 


Component 

A002 

A201 

A203 

A212 

A210 

A410 

A  602 

A803 

A601 

A810 

A800 

A603 

A400 

A613 

A  802 

A801 

A412 

A4 14 

A413 

A204 

A411 

A412 

A012 

A612 

A011 

A611 

A214 

AO  14 

P203 

P210 

P310 


Va  lue 

0.6000 
-0.2500 
-0.3500 
0.3890 
0.0500 
-0.9999 
-0.5000 
-0.3000 
0.2500 
0.1000 
-0.9999 
-0.5000 
0.2500 
-0.5000 
-0.2500 
0.2500 
0.0500 
0.0500 
C  9000 
C. 3500 
-0.5439 
-0.0500 
-0.6244 
-0.4381 
0.2500 
0.3535 
-0.0000 
-0.4000 
-0.0389 
0.0700 
-0.0500 
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STATIC  CHECK 


Component 

P601 

P810 

P701 

P9I0 

P603 

P513 

P612 

P012 

P711 

P311 

P213 

P311 

P312 

T311 

T312 

T303 

T310 

T300 

T301 

T313 

T400 

T302 


Va  lue 

-0.0144 
-0.0043 
-0  2500 
-0.1000 
-0.1000 
40.0244 
40.0381 
-0.0624 
-0.0438 
40.4000 
40.0010 
'0.0050 
40.0100 
40.2500 
-0.6000 
40.5000 
40.3000 
40.5000 
40.6000 
40.3486 
40.5439 
40.4000 
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SWITCH  ASSIGNMENT 
(AUTOMATIC) 


Sense  Line 

Exc 

Function 

C 

SW  1012 

Print  Out  Check  &  Time  Hist 

*1 

SW  1211(1) 

P£  *  f  (Pr,  P0,  u) 

*2 

SW  1053 

(0) 

P0  =  f(Pr,  P§,  (i) 

3 

SW  1013 

Not  Used 

4 

Used  as 

Monitor 

(for  H<-10km> 

5 

Used  as 

Monitor 

(for  H^lOkm:  Scale  Change) 

*6 

OP 

Hold  Simulated 

*7 

IC 

IC  Simulated 

*AstcrJsks  denote  sense  lines  which  are  hard  wired 


I IA3 .  SYSTEM  FLOW  DIAGRAM 


Figure  (II A3 j  represents  the  flowchart  of  the  system.  The  purpose  of 
this  diagram  is  to  shew  the  flow  of  the  convergence  scheme.  For 
derivation  of  and  insight  into  the  ray  tracing  equations  see  AFCRL 
Report  A69-0123 ,  Contract  No.  FI 9  628 -67-C0359,  February  1969. 

Included  also  is  the  digital  listing  of  the  program. 


Figure  (ITA3.)  System  Flowchart,  Cent. 


73 


SUBROUTINE  ANALOG  SE1  -UP 


SUBROUTINE  RUN  DATA 


Figure  (IIA4.)  Subroutines  Analog  Sst-Up  and  Run  Data 
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'  -  mm 


Figure  (IIA5.)  Subroutine  Initial 


I 


D*s[2(l-X)-Y*Sln*^±[(YSin  ^)4  +4Yt(l-X)*Co*#^]^  ]* 


tftSln\^^l-2X(l-xyin^ 


1  |£.|2XH-X(l-X)j2+4(l-X)tYCo«y)  \ ^/|D^I-2X(I-X) j. 

l5^s|x(l-X)Y|’-2(l-Cos  V)  +  i(2(YSin^)*(l-CosV)+4(l-X),Cc*  2X(I— X) 

i  |^s  j2X(l-X)Y*Sinf  Co»*  pYSIn^)*-  2(l-X)*-lj  Jy^l-  2X0-X^j 


3*  PrCo«*-7*Nr 


3ps 

3* 

P^Coa  i 

3 PS  M  Sin  f 


l  3m  I  3u  3N 


RETURN) 


Figure  (IIA7.)  Subroutine  Partial  MU 

REF:  fPP  3-24,  3-25  AFCRL  69-0123) 


P 


■ft 


START 


h-hinin 


h^mo* 


Figure  (IIA8.)  Subroutine  Electron  Density 

REF.  (Pg  (3-22)  AFCRL  69-0123) 


tir> 


START 


Figure  (IIA10.)  Subroutine  Set  Pot 


A=Cosift$  SlrSfl 
B'Sin^gSinflg 


Cosf>M=  O74022A-.I86I278B+.97992  Co*  8q 
Sin  0M  =  -/l-Cos20M 


>  /  \  < 

- CTcosSm'O^ — 


TAN"1  (  It  37.29578 

V  Cos  8ft' 


STM  =TAW-1  (  37.29578+180 


Sin  <f>M:  (.93358A  +  .33837B ) 

Sin  flu 

Cos  (.33tl739A— .914833 78 -.19937 Cos  8Q ) 

Sin«?M 


>  /  ,  X  < 

— — — <x.  Cos  0  _p - 


— *  $>M3  bt *  380 

\ .  - , - 


Figure  (IIAli.)  Subroutine  GEOTQH 

REF:  (P3-21,  AKCRL  69-0123) 


IIA3.2  DIGITAL  PROGRAM  LISTING 


SFO». TUI , *4021 

DIMENSION  ANGI  i C (  3 ) ,R I NPUT ( 8 ) .OUTPUT ( 13 ! *SF(8) 

DIMENSION  OUT! ( 100) ,0UT2 ( 100 ) .OUT 3 ( 1 00 ) ,0UT4( 100) 

DIMENSION  0UT5 (  100  ) 

DIMENSION  ALI3) *  THF  <  6 ) ,HERR(6) 

COMMON /SCOOP /HGT ( 100) , ED ( 1 00 ) , KMAX ,HMAX , HTOP 

C0MM0N/SC00P1/SCALI7) ,P0TADR(3I 

COMMON/ SCOOP? /H  , THO ,PH I  0 .ALPHA , BETA »  IRON 

COMMON/ SCOOP 3/ A 1 1.A12.A13.A21 ,A72,A23*A31 »A32*A33 

COMMON /SC00P4 /PRO ,PTH0,PP0 

COMMON/ I NFO/FKCTtFKC* SIGN 

common/ i nfoi /r.costh, si  nth, cospm.sinpm.thd.ph id, opt, sop t 

COMMON/ ! NF02/Y  »  YR *  Y  TH » YSQ.NR  »MTH *X  »XR »OMX  »OMX 2  » XOMX 
COMMON/ I NF03/MUCK , PEN , PTHN , PPN ,COSPS *C0SP2  » S INP2 ,S I NPS 
COMMON/ INF04/YL * YT , VL2 , V T 2  * MU.MUSQ, Z  * Z2  * 2R *ZMOD 

COMMON/  INF0)5/PSpR»PSPTH«PSPP»P.4I  TH*MUDMU  »MURMU  ,MIJTMH,MUPHM»OPMUPO 
COM''  )N  /  I NFOS  /RM  ♦M?*AC*A1»A2»A4«AS  ,  A6  ,  A7  ♦  Aft  ,  BO  ,9 1 ,02 , 94  ,R?  »B6 , !37 ,  R8 
COMMON/ INF07/D3  »D2  >  A  *B  »R S4 ,RS ,THS  ,S1 ,  '  2 , DSQ  *  AA  » AA2 , 0B ,092 
COMMON/  INFOS /R'« 4* ThRM, Ml  »MUZMU  .MLJXMU  *MUVMU 
COMMON /COORD/THTG  , PHG  *  THTM , PHM 
EXTENDED  POTADR 

REAL  NR  *MTH  »Ml)CK  »MIJ  *MUSQ  »MUPMU»MURM()  »MUTMU  ,MUPHM  *M2  *  MUZ  MU » MUX  MU* 
1MUYMU »N 

LOGICAL  INDICL. INDICM, INF ICN, INDICK 
LOGICAL  INTFG8 

LOGICAL  INTEG*tNTEGI » INTFG2  » I DAVE 
LOGICAL  I H  1ST  t PN'>, NOTBE 

EQUIVALENCE  (SNPT] ,H)  , (SNPT2.TH) » (SNPT3*PHI ) 

EQUIVALENCE  iSNDT4,PR) , ( SNPT 5  * PTH) *!SNPT6*PP) 

EQUIVALENCE  (SNPT7.P) 

EQUIVALENCE  (SNPT8.D) 

DATA  AMPL/4HAOOO/ 

DATA  A1I.A12.A13/  0.3M1  739,-0.9148337,-0.19937/ 

DATA  A? 1  *  A??  » A2  */  0 . 93 388 , 0 . 0 58 37 , 0. 0/ 

DATA  A31.A32.A33/  0.0714322,-0.1861278,0.97992/ 

DATA  (SF(J).J=1,8) /loop. 0,3. 4889 *6. 9770 *2. 0.2.0, 2.0*5000,0. 

6?  DO').  0/ 

/  LDOR  =*40040, ,  1? 

CSJ.$S«^SSDATA  READ  IN  AND  VF:R  I F  I C  A  T  I  ON  PRINT  OUT  OF  ELECTRON  DENSITY 
CALL  DATA  READ  IN 

C*»*#**#*REaO  SCALE  FACTORS,  °0T  ADDRESSES 
CALL  analog  SETUP 

C##*#*##*Rf AD  IN  RUN  DATA  AND  OR  I  NT  OUT 
8400  CALL  RUN  DATA 
C***#»*#*RESET  SFNSE  LINFS  0,1, ?,3 
no  5  6  K  =  1 f 4 

j  =  K-l 

CALL  OTE^FO? 1 ,J.  INOICL, IERR) 

88  '*0N7INUr 

C*#**#**»  INITIALIZE  STATE  VARIABLES  THRU  COORDINATE  XFORMAT  ION 
C*#INIT IALIZE  TIME  HI  STOP  t  ARRAYS 

I  COUNT *0.0 

E  NO* ♦ F AL  SF • 

A'OTBF  =  .TRUE  . 


ftHlN*  ms  Mut 
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TAU=0.0 
I TEST=1 
STFP=10.0 
ILAST=100 
0UT1 ( 1  )  =0.0 
OUT2 ( 1 ) =0.0 
OUT3 ( 1 ) =0 . 0 
OUT  4 ( 1 ) =0.0 

C**0EGIN  ALPHA  APPROX  CALCULAT I ON 

READ  (6,610)  THTARG. PUT ARG .DTARG, THTRAN .PHTRAN .DALPH. ORE TA» CONST 
WRITE  (6.610)  THT ARG, PHT ARG. DTARG, THTRAN .PHTR AN .DALPH.OBET A .CONS! 

610  FORMAT t 6X .2F10.4 ,F 5 . 1 .2F10. 4 .2F5. 1 .F 10 . 3  ) 

8ETA=221 .43 
BOR IG=BFT  A 
BE TAE= 2 7 0.0-BET A 
TOL 1=0. 01 
TOL2=10.O 
ZMOD=0.0 
1=0 

616  1=1+1 

7979  FORMAT ( 6H  I »  5X (5HALPHA.5X.5H  BETA.5X.5H  H  .5X.5H  THD  .5X. 

15H  PHID//  ) 

I F ( I.GT.3)  GO  TO  700 
ALPHA=12.5*FLOAT(  I  ) 

7777  FORMAT ( 5X  » F8  •  2  ) 

AL ( I ) =ALPHA 

660  THO=THTRAN 
PHIOsPHTRAN 
CALL  INITIAL 
I T IME=0 
HZERO=0 • 0 

C****#**#CALCULATE  AND  SET  IC  POTS.H » THETAM  »PH IM 
POT1=HZF.RO/1000.0 
POT2=THTM/200.0 
POT3=PHM/400.0 
CALL  QPSO(l-.IERR) 

CALL  SETPOT  ( POTADR ( 1 ) , POT  1 ) 

CALL  SETPOT  ( POTADR ( 2 } ,POT2  ) 

CALL  SETPOT ( POT  ADR ( 3 ) » POT  3 ) 

CALL  ORDALO ( 1 »AMPL .DUMMY  » I  ERR ) 

C#**#»#**CALCULATE  AND  TRANSFER  IC'S  PR.PTH.PP 
PR=PR0 
PTH=PTH0 
PP=PP0 
TH=TH0 
PH  I =PH I  0 
TMO=THO 
PM0«PHI0 
TM»TMO 
PM=PMO 
S«0.0 
TAU«0.0 

C#*#*#***CALCULATE  INDEX  of  REFRACTION 

CALL  INDEX  OF  REFRACT  ION ( H. TH. PHI .PR.PTH.PP ) 
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PRO=PRO*MU 
PTHO=PTHO*MU 
PP0  =  PP0*MLJ 

C**#*#**#SET  UP  BUFFER  ARRAY  FOR  ICS 
ANGLICm  =-PR0/2.0 
ANGLICf?) s-PTHO/2.0 
ANGLIC  3)  =-PP0/2.0 
C*#*###**LOAD  AND  TRANSFER  DAC'S 
2228  CALL  ODALDO ( 1 3 » 3 » ANGL IC ♦ IERR , IECHL ) 

IF ( IERR.FQ.2 )  WRITE<6,81)  IECHL 
IF{ IERR.FO.3i  WR I TE ( 6 , 82  ) 

CALL  ODAXRO(12»4,IERR) 

I F  C IERR.EQ.2)  WR ITF ( 6  ,  S3  ) 

CALL  QICOU  .IERR) 

CALL  CTFFF0(1*7  .INTEG2.1EE) 

CALL  OSYNCO (1.7, INTEG2  ) 

I  PA SSI  =  1 
GO  TO  70 
11)8  IPASS1=0 

WRITF(6»772S)ALPHA,BETAE 

TYPE  1108 

PAUSE 

C*****T  EST  SLO  FOR  PRINTOUT 

CALL  OTEFFO; 1,0, I  DAVE, IERR) 

IF  I  I  DAVE)  GO  TO  70 
1936  I  DAVE=. FALSE . 

TYPE  1109 

11.08  FORMAT  (  42H  BLIP  F$W  101?  FOR  If  PR  I  NT ...  PRESS  FLAG  8) 
.1100  FORMAT  ( 26H  PRESS  FLAG  8  TO  CONTINUE.) 

PAUSE 

/  SFL  = '  63 ,  »  0 

/  LOT  =>77777 

CALL  0000(1. IERR) 

CALL  QT  FFFO ( 1 .6  » INTFG1 . IFE  > 

CALL  OSYNCO ( 1 .6. INTFG1  ) 

C****** **ENTER  HYBRID  LOOP 
fCi  CONTINUE 
/  STT  / 1  TIME 

/  LOT  ='77777 

SFL  =  '  6  2  ♦ » 0 

6662  DEL  TAT- (32767-1  TIME) *SCAL( 7) 

(;#**#*#**RE4r>  ADC'S 

CALL  OSTRF0(O.16.IPRR) 

T F { IFRR.EQ.2)  WRITE  <6,00) 

CALL  OA OF VO ( 0  >8  »  R ! NPUT . IERR  ,  IECHAN  ) 

I F ( IERR.FQ.2)  VRITF(6,91)  IECHAN 
!F(  IERR.EQ.3)  UR ! TE ( 6 . 92  ) 

CALL  OTRCKO  (0,16. IERR  ) 

I F ( IERR.EQ.2)  WRITE  (6,93) 

r*###***#ALL  A N A  LOG  INPUTS  RECEIVED  -  UNSCALE  ANALOG  INPUTS 
C*#«CI!FCK  FOR  H  EQUAL  1C.  IF  SO  CHANGE  SF(1) 

CALL  OT  F  c  F  D ( 1 ♦ S  » T  N  T  E  C  8  » IFE) 

CALL  OTEFFO ( 1 .5 , IN  TEGS  » IERR ) 

I F { INTFGS ) SF ( 1 ) =19.0 


C*#***CHECK  FO  H  LT  -10  KM 

CALL  QTEFF0(1»4»INTEG,IEE) 

I  F (  INTFG)SF(  1)=1000.0 
SNPT1=SF( 1 )*RINPUTt 1) 

SNPT2=SC( 2 ) *R I NPUT ( 2) 

SNPT3=SF ( 3) #R INPUT ( 3 ) 

5NPT4=SF ( 4 ) *RI NPUT ( 4 ) 

SNPT5=SF ( 3 ) *R I NPUT ( 3 ) 

SNPT6=SF(6)*RINPUT (6) 

SNPT7=SF  ( 7 ) #R I NPUT  < 7 ) 

SNPT8=SF ( 8 ) *R I NPUT ( 8 ) 

SFtl) *1000.0 
30  CONTINUE 

C****#*#*»  TEST  SL6  FOR  OPERATE  MODF  OF  ANALOG*  IF  ANALOG  NOT  IN  THE 
C#**##****  OPERATE  MODE  GET  DELTAT*0.0 
CALL  0TEFF0(1»6*INTEG1»IEE) 

IF ( .NOT • INTEG1 )  DEL  TAT  =  0.0 

CERC  =  SQRT(DTHC*DTHC+SINTH*SINTH*(DPC*DPC  )  ) 

tm=tm+deltat*dthc 

PM=PM+DELTAT*DPC 

TH  =  TM 

PHI=PM 

S=S+6370.0*CERC*DELTAT 

C********»  TEST  SL7  FOT  IC  MODE  OF  ANALOG*  IF  ANALOG  IN  IC  REINITIALIZE 
C*********  THETA  AND  PHI 

CALL  QTEFF0<1»7,:nTEG2,IEE) 

IF( .NOT. INTEG2 )  GO  TO  8763 
TM=TMO 
PM=PMO 
3  =  0.0 
TAU=0.0 
8763  CONTINUE 

IF ( .NOT .ENDJGO  TO  71 
T  AU=TAU+DELTAT 
I  COUNT *T  AIJ/STEP +0.001 
IF ( ICOUNT .GT.ILASTJGO  TO  71 
I F ( ( ICOUNT+1 ) .EO. ITEST )GO  TO  71 
ICP  * ICOUNT+1 
OUT  It ICP) =TAU 
OUT  2 ( ICP) =H 
OUT3( ICP) *S 
01JT4(  ICP)  *D 
OUT5( ICP) ~P 
ITEST*ICP 
71  CONTINUE 

C»*##**#*CAl CUL ATE  MUtELECTRON  DENS  I TY t ELEC TRON  DENSITY  DERIVATIVES* 
C######*#MAGNFTIC  FIELD, AND  MAGNETIC  FIELD  DERIVATIVES 
CAuL  INDEX  OF  REFRACT  I  ON ( H* TH.PHI »PR »PTH *PP ) 

C***##***APPL I CAT I ON  OF  CONSTRAINT  BETWEEN  MU  AND  PRtPTHfPP 
C*##*#**#FOR  DETERMINING  CORRECT  PR*PTH,PP 
C##***#*#UNNORMAlIZE  PR 
PR«PRN*MU 

C####>mm**TEST  SL1  TO  CHOOSE  PP  AS  A  FUNCTION  OF  PR, PTH. MU 
CALL  CTEFF0U.1  .  INDICL  .KERRSS) 
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IF  (KERRSS.EQ.2 )  PRINT  4001 
IF  (KERRSS.EG.3 )  PRINT  4002 
IF ( .NOT. INDICL )  GO  TO  969 
C*#*##***HERE  WE  HAVE  CHOSEN  PP  =  F ( PR  *  PTH *MU ) 

PSIGN=-1.0 

IF  (PP.GT.O.O)  PSIGN=1.0 
RADCL*MUSQ*( l.-°THN«PTHN)-PR*PR 
ARADCL=A.GS(RADCL) 

PP=PSIGN*  SORT(ARAOCL) 

GO  TO  1Q69 

C****##**TEST  SL2  TO  CHOOSE  PTH  AS  A  FUNCTION  OF  PR»PP*MU 
969  CALL  QTEFF0(1.2tINDICM,ICRUS> 

I F ( ICRUS.EQ.2)  PRINT  4001 
IF! ICRUS.EQ.3)  PRINT  4002 
IF(.NOT.INOICM)  GO  TO  1969 
C»##**#*#HERF  WE  HAVE  CHOSEN  PTH  =  F  fl-' K  *PP *MU  ) 

RPSIGN=-1.0 

IF(PTH.GT.O.O)  RPSIGN=1.0 
RADCL'I  =MUSQ*  (1«-PPN*PPN)—PR*PR 
ARADIC=ABS(RADCLI) 

PTH=RP5IGN*SQRT ( ARADIC ) 

1969  CONTINUF 

C*#*###**PR,PTH*PP  ARE  NOW  DETERMINED*  WITH  IR  WITHOUT  THE 
C###**#*# APPLICATION  OF  THE  CONSTRAINT 

C***-:<#*##CALCULATION  OF  PARTIAL  DERIVATIVES  OF  MU  NOW  FOLLOW 
CALL  PARTIAL  MU ( PR  ,  PTH *PP ) 

C****#***CALCULATION  OF  STATE  VARIABLES  FOLLOW 
C*#******THETA  DOT/C 

DTHC= (FTH/MUS0-MUPMU#PSPTH> /R 
C********PHI  DOT/C 

OPC=(PP/MU50“MUPMU*PSPP) /( R*$INTH) 

C#*#*#***PR  DOT/C 

DPRC  =  MtJRMU  +  PTH*DTHC  +  PP*  DPC»SINTH 
C********PTH  DOT/C  +  PTH*R  DOT/C 

PTHOTM=MUTMU/R  +  PP*COSTH*DPC 
C********PPHI  DOT/C  +  PPH I *R  DOT/RC 

PPHOTM= (MUPHM/R  -  PP*COSTH*DTHC )  /  SlNTH 
C* * **•*■**# LOAD  DAC  BUFFERS 
GR01JP  =  1  ,+OPMUPO/MU 
OUTPUT ( 1 ) =SCAL  ! 1 ) *5O.0»DPRC 
OUTPUT ( 2 ) =SCAL <  2) *500.0*PTHOTM 
OUTPUT (3 ) *SCAL (  3 ) *  500 . 0*PPHDTM 
OUTPUT' 4) =SCAL (4)*DTHC*1000.0 
OUTPUT ( 5 ) =SCAL ( 5 )*DPC* 1000.0 
OUTPUT! 6) = (PR/MUSQ) /2.0 
OUTPUT ( 7'*-SCAL (6)*MUPMU*PSPR/2.0 
OUTPUT ( H ) =SOOO. 0»PTH/R 
OUTPUT!  9)  =  V)00.0*PP/R 
OUTPUT! 10 ) s-MU/2.0 
OUTPUT! 1 1 ) =DELT AT/ 100.0 
OUTPUT! 12 ) =S/1000.0 
OUTPUT! 1 3 i »-GROUP/ 50. 0 
C*###*#*#TEST  SLO  FOR  PRINT  OUT  CHECK 
CALL  OTEFFO!1.0,INOICK*IERSS) 
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IF{ IERSS.EQ.2)  PRINT  400] 

IF  C IERSS.EQ.3)  PRINT  4002 

IF(. NOT. ( INDICK. OR. IDAVE > )  GO  TO  506 

C*#****##IF  INDICK  TRUE  (SLO  SET)  WE  PRINT  OUT 
N*X /FKCT 
DNR=XR/FKCT 

C#**##**#PRINTING  FORMATS 
WR  I  TE ( 6 .77  ) 

WRITE  (6*9140)  PR,PTH,PP 
WRITE  (6*9141)  H.THO.PHID 
WRITE  (6*914/)  NR.MTH.COSPS 
WR I TE(6 ,9145 )  Y.YR.X 
WRITE  (6*9153)  XR.Dl.D? 

WRITE  (6*9700)  Z.2R.Z2 
WRITE  (6*9701)  A.B.RS4 
WRITE  (6*9702)  RS.THS.Sl 
WRITE  (6,9703)  S2,DS0*AA 
WRITE  (6,9704)  AA2*BB,BB2 
WRITE  (6*9705)  RM4  »RM»THRM 
WRITE  (6*9706)  M1,M2»MUSQ 
WRITE  (6*9707)  A0,A1,A2 
WRITE  (6*9708*  A4,A5,A6 
WRITE  (6*9709)  A7»A8,0O 

WRITE  (6,9710)  B1»B?,B4 
WRITE  (6*9711)  B5  *  R6  *  07 
WRITE  (6*9712)  B8 »  PS  I  TH  ,MUXMtJ 
WRITE  (6*9713)  MUYMU»MUPMU*MUZMlJ 
WRITE (6, 91 59)  PTHOTM.PPHDTM .GROUP 
WRITE  (6,9156)  PSPR *P$PTH , PSPP 
WRITE  (6*9157)  MlJRMU  *MUTMU » MUPHM 
WRITE  (6,9158)  DTHC.DPC.DPRC 
WRITE  (6*9154)  MU, MUCK 
WRITE(6,8150)  COSTH.SINTH 
WRITE (6*8151)  OPT.SOPT 
WR I TE ( 6  *8 1 52  )  N.DNR 
WRITE (6*9155)  TH.PHI 
WR I TE ( 6  *77  ) 

WRITE(6,916) 

916  FORMAT ( 37H  J  OUTPUT(J)  CHANNEL//) 

DO  816  J=l*13 
JCH* J  »l 

WRITE  (6*9160)  J .OUTPUT ( J )» JCH 

9160  FORMAT( T3,6X*F14.6, 1 1X.I2  ) 

816  CONTINUE 

WR I TE ( 6 , 77 ) 

WR I TE ( 6*9922 ) 

9922  FORMAT ( 37H  J  RINPUT(J)  CHANNEL//) 

DO  47  J  =  i  ,6 
JCH-J-l 


WR I TE ( 6  *9160 )  J ,RI NPUT ( J ) * JCH 
47  CONTINUE 
PAUSE 

C**##*##»DAC  LOAD  AND  TRANSFER 

506  CALL  QOALDOIO. 13, OUTPUT, IERR4.IECHU 


IF( IERR4.EQ.2)  PRINT  81.IECHL 
I  F (  IERR4.EQ.3)  PRINT  82 
CALL  QDAXR0(0*16,IFRR5) 

I F ( IERR5.EQ.2)  PRINT  83 
C#*******TEST  SL3  FOR  NEW  RUN 

CALL  QTEFFOI 1 *3  * INDICN, IERR) 

I F (  INDICN)  GO  TO  8400 
I F ( I  DAVE ) GO  TO  1936 
D  I  EF  =  ABS ( PHTARG-PH I D ) 

IF  C  I  PASS  1  •  EQ.'l  )  GO  TO  1118 
IFIDIFF.GT.TOLl )  GO  TO  70 
CALL  QHOLDO ( 1 ♦ I  ERR ) 

I F ( END )GO  TO  913 
THE (  I  ) =THD 
TH0RN=90.0-THD 
HERR ( I ) = 1 0. 0*H 
I F (  I.LE.3)  GO  TO  7982 

7982  WR I TE ( 6  *  7979 ) 

WRITE (6. 798 3) I , ALPHA , BETAE ,H .THORN , PHI D 

7983  FORMAT ( 1 5 , ?X ,F8 . 2 , ?X . F8 . 2 . 2X .F8 . 2 .2X .F8 . 2 .2X »F8. 2/ ) 
f'NFW=HERR(  I  ) +ABS ( CONST# ( THE (  D-THTARG)  ) 

D  I  FFD  =  ABS ( DNFW ) 

IF (DIFFD.LE.TOL2)  GO  TO  6b3 
I F ( I -3 )  615.690.700 
690  HOR IG=20000.0 
DO  695  K=  1 .3 

IF (HERR (K ) .LT.HORIG)  GO  TO  692 
GO  TO  695 
69?  AORIG=AL(K) 

HOR I G=HFRR ( K ) 

THOR  IG=THE  (  K  ) 

695  CONTINUE 
651  1=4 

ALPHA=AOR I G+DALPH 
BET  A  =  ROR I G 
GO  TO  650 

70n  I F ( 1-5)1701 ,1703.1703 

1701  DORIG=HORIG+AR5(CONST*{THORIG-THTARG) ) 

DADA  =  HERR(  I  ) +CONST* ( THE (  I ) -THTARG ) 

PDPA= (DADA-DORIG) /DALPH 

C*»*  TEST  LOGICAL  VARIARLE  NOTBE  TO  DETERMINE  8ETA  OMISSION 
I K ( .NOT. NOTRE )  GO  TO  1754 
PDPRF  =  0.0 
GO  TO  1753 
1754  1=5 

RFTAsRORJG+DRETA 

alpha=aorig 

GO  TO  650 

1702  DBDB  =  HERR (  I ) +CONST * ( THE ( I ) -THT ARG) 

POPR5MDBDB-DORIG)  /DBE7 A 

1763  GRADD?=PDPA*PDPA+PDPBE*PDPBE 
GRADD-SQRT (GRADD?) 
dalphi =DALPH 
adalfi^arsidalphi ) 


DALPH=-(DADA/GRADD2)*PDPA 
ADALF=ABS ( OALPH ) 

KAD=  1 

IFCAOALF.GT.ADALFl )  KAD=APALF/ADALF1 
DALPH=DALPH/KAD 
D8ETA=-  (DADA/GRADD2 ) *PDPBE 
C*#*##*#*  TFST  NEW  ALPHA. BETA 
ALPHA=ALPHA+DALPH 
RETA=BORIG+DBET  A 
HORlG=HERR< I  ) 

1  =  1  +  1 
GO  TO  6*50 

C*«**«*«*  RF-INITIALIZE  FOR  CONTINUED  SEARCH 
1703  aorig=alpha 
BOR  I G=BET A 
THORIG=THD 
I F ( END ) GO  TO  650 
IF(NOTBE)  GO  TO  1701 
GO  TO  651 

C********  STATEMENT  725  INDICATES  SUCCESSFUL  TARGETING 
663  TYPE  725 
FND=. TRUE  • 

725  FORMAT ( 2 9H  THIS  IS  THE  OPTIMUM  RAY  PATH) 

7725  FORMAT ( 7H  ALPHA= ,F8 . 2 , 1  OX . 5HRET A= . F8. 2/ ) 

PAUSE 

GO  TO  1703 

C**SET  SLO  FOR  TIME  HISTORY  PRINTOUT 

913  TYPE  916 

915  FORMAT ( 25H  SET  SLO  FOR  TIME  HISTORY) 

PAUSE 

C*#*T IME  HISTORY  ARRAY  CHECK  FOR  PRINTOUT 
CALL  QTEFFOd.O.IHIST.IERR) 

I F ( IH I  ST ) GO  TO  914 
926  FND=. FALSE. 

CALL  QICO(l.IERR) 

GO  TO  70 

C*###PR  I NTOIJ  i  OF  TIME  HISTORY 

914  WRITE' 6.2  121) 

2121  FORMAT! 10X.15H  H  =  ALT I TUDE ( KM ) / ♦ 10X , 28H  S=GREAT  CIRCLE  D I  STANCE (KM ) 
6/.10X.17H  D=PHASE  PATH! KM ) / v 10X . 1 7H  P=GROUP  PATH(KM)/) 

DO  89  ICP=1 » I  TEST 

WRITE (6, 80) I  CP. OUT  1 ( I  CP) *OUT2( I  CP) 

WRITE (6. 84,  OUT  3! I  CP) .0UT4! I  CP) .OUT5 ( I  CP) 

80  FORMAT ( 5H  iCP*I2.5H  T AU=F 10. 2 ♦ 10X . 3H  H*F10.6) 

84  FORMAT { RX  «3H  S  =  F  10 . 6 . 10X . 3H  D  =  F 10 . 6 . 1  OX ♦ 3H  P  =  F10.6) 

89  CONTINUE 
TYPE  444 

444  FORMAT ( 32H  TIME  HISTORY  PRINTED.  RESET  SLO) 

PAUSE 
GO  TO  926 

C *»#### ##FORM AT  STATEMENTS 
90  FORMAT  (31H  BLK  ADR  ERROR  !N  STORE  ROUTINE) 

91  FORMAT  (23H  ADC  CHANNEL  OVERLOAD  -.15) 

92  FORMAT  (21H  NOW-EXISTING  CHANNEL) 
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93  FORMAT  H1H  BLK  ADR  ERROR  IN  TRACK  ROUTINE) 

81  FORMAT  (25H  OVERFLOW  IN  DAC  CHANNEL  *15) 

82  FORMAT  (25H  NON-EXISTING  DAC  CHANNEL) 

83  FORMAT  (22H  DAXFR  BLOCK  ADR  ERROR) 

'001  FORMAT ( 2 1H  CONSOLE  DISCONNECTED) 

4002  FORMAT ( 33H  NON-EXISTING  COMPONENT  REQUESTED) 

77  FORMAT! 1 H 1 ) 

8150  FORMAT ( 7H  COSTH= »F 7 . A , 1 1 X « 7H  SINTH=*F7.4) 

8151  FORMAT ( 5H  OP T = * E8 . 4 . 1 2 X , 6H  S0PT=*F8.4) 

8152  FORMAT ( 3H  N = * E 1 3 . A  * 9X  ,  5H  DNR=*E13.6) 

9155  E ORMAT ( 4H  TH= * F 10 . 6  ♦  1 1 X  * 5H  PHI**F10.6) 

9147  FORMAT ( 4H  NR = * F 7 . 5 , 1 4X * 5H  MTH= . F7 . 5  *  1 3X * 7H  COSPS=»F7.5) 

9140  FORMAT ( AH  PR= , F 7 . 5  *  1  AX , 5H  PTH= » F7 . 5 ♦ 1 3X *  AH  PP=*F7.5) 

9141  FORMAT ( 3H  H= * F 1 2 . A  *  1 OX ♦ 5H  THD- »F12.6*8X*6H  PHID=»F12.6> 

9145  FORMAT ( 3H  Y=*F7»5*15X*AH  YR =  .  F9 . 5 . 1 2X , 3H  X=*F10.5) 

9153  FORMAT I 4H  XR  =  , E 1 1 . 6  * 8X , 4H  D 1  =  * E 1 3 . 6  *  8 X , AH  D2=»E13.6) 

9154  F  ORMAT  I  4H  MU=  *  F  7. 5 , 1  AX  , &H  MIJCK=*F7.5) 

9156  FORMAT ( AH  PSPR= *F 1 2 .6 . 7X ♦ 7H  PSPTH= *F 1 2. 6 *6X »6H  PSPP=»F12.6) 

9157  FORMAT  ( 7H  MURMU= , F 1 2 . 6 . 6X . 7H  MUTMU= . F 1 2 .6 » 6X ♦ 7H  MUPHM= *F12. 6 ) 

9158  FORMATIAH  DTHC=  . F 1 2 . 6 . 7X * 5H  DPC  =»F12*6*8X»6H  DPRC=*F12.6) 

9159  FORMAT  (  8H  P  THOT  M  =  ,  F 1  2 . 6 . 5X  » 8H  PPHDTM.=  *F  1 2 .6  » 5X  ,  7H  GROUP=*F7.4) 

9700  FORMAT ( 3H  Z = . E 1 3 . A » 9X , AH  ZR = * E 1 3 . 6 » 8X * AH  Z2=*E13.6) 

9701  FORMAT ( 3  H  A = , E 1 3 . A  * 9X  .  3H  R = , E 1 3 .6 . 9X  ♦  5H  RS4=*E13.6) 

970?  FORMAT  4H  R r = , E 1 8 . 6 ♦ 8X , 5H  THS* * E 1 3 . A . 7X . AH  S1*.E13.6) 

9703  FORMAT (4H  S2= *E  13. 6 *8X  ,5H  DSO= ♦ E 1 3.6 . 7X ♦ AH  AA=,E13.6) 

9704  FORMAT ( 5H  AA2= » E 13 .6 . 7X ,AH  BB= , E 1 3 . 6 . 8X  .  5H  BB2=.E13.6) 

9705  FORMAT ( 5H  RMA= , E 1 3 . 6 , 7 X , AH  RM= ♦ E 1 3. 6 » 8X ♦ 6H  THRM=,E13.6) 

9706  FORMAT I AH  Ml = . E 1 3. 6 » 8X , AH  M2= * E 1 3 . 6 . 8X , 6H  MUSU=*E13.6) 

9707  FORMAT ( 4H  A 0= , E 13 . 6  * 8X , AH  A  1  =  » F 1 3 . 6  * 8X . AH  A2=»E13.6) 

9704  FORMAT ( 4H  A4=,F13.A*8X,AH  A 5 r * E 1 3 . 6 , 8X , AH  AA=*E13.6) 

0709  FORMAT ( 4H  A 7 = , E 1 3 . A . 3X , AH  A 8 - » E 1 3 . 6  * 8X » AH  B0=*E13.6) 

9710  FORMATIAH  R1=,E13.6*8X,AH  B2* » E 1 3 . 6 ♦ 8X , 4H  BA=»E13.6) 

9/11  FORMATIAH  R5= *E 13.6 *8X *4H  B6  =  * E 1 3 . 6  * 8X » 4H  B7=»E13,6) 

9712  FORMATIAH  R8 = , E 1 3. A . 8X » 7H  PS  I TH= »E 1 3 . 6  * 5X *7K  MUXMU* . E 1 3 . 6 ) 

9713  FORMAT  (  7H  MU YMU=  * E  1 3 .6 . 5X  .  7H  MIJPMU*  »E  1 3 .6  » 5X  »  7H  MUZMU*  .E13.6  ) 
END 
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3158  FORMAT (6H  DTHC  = , FI  2. 6, 7X ,5H  0PC= ,F12 . 6 , 8X, 6H  DPRC=,F12.6> 

9159  FORMAT OH  PTHOTM= , F12. 6, 5X ,8H  PPHOTM=, FI  2. 6, 5X , 7H  GROUP  =  ,F7.4> 
ENO 


SUBROUTI ME  ANALOG  SETUP 
COMMON/SCOOP  1/ SC AL ( 7 ) , POTA OR  (3 ) 
EXTENDED  POT  ADR 
READ (5 ,77)  <  SC  AL  OO  ,  K=  1,  7) 

NRI T  £  <  6i  77MSCAL(K),K-=1,7> 

77  FORMAT (7F11. 5) 

READ  (9 ,7/7  >  (POTAOR(J) ,J=1,3I 
777  FORMAT  (A4,4X, A4.4X.A 4) 

RETURN 

END 


SUBROUTINE  RIJN  DATA 
COMMON/INFO/FKCT,FKC,SIGN 
COMMON/SCOOP2/H, THO , PH  10 , ALPHA , BETA , IRUN 
8400  RE  AO (5 , 3  4)  IRUN 
34  FORMAT (18) 

RE  AD (5, 7033) H, TMO , PM  10 , ALPHA , BET  A 
7033  FORMAT  (5F15. 7) 

REAO(S,7033)  SIGN, FKC 
WRITE (6, 7000) 

7000  FORMAT (1H1) 

WRITE(6, 7010)  IRUN 

7010  FORMAT  T12H  RUN  NUMBER  ,13//) 

HRITECG./Oll)  SIGN , FKC 

7011  FORMAT (6H  SIGN=, F4. 1 ,6 X, 7HFREQKC=, F8 . 2/) 
HRITE(6,7012)  H,THO,PHIO 

70120FORMATf3H  H= ,F 7 . 2 , 6X , 7HTHETAG* , Fb. 2, 3X , 5HPHIG* , 
1F6.2/) 

WRITE! 6, 7013)  ALPHA, BETA 
7013  FORMAT (7 H  AL PHA= ,Ffe . 2, 3X , 5HHET A- ,F6. 2/ > 

RETURN 

END 


C 


SUBROUTINE  INITIAL 

COMMON/INFO/FKCT , FKC, SIGN 

COMMON/SCOOP  2/ H, THO, PH  10, ALPHA, BET  A, IRUN 

COMMON/SCOOP  3/All,Al2,A13,A21,A22,A23,A31,A32,A33 

COMM ON /SCO OP 4/ PRO , PT HO ,PP0 

COMMON/COORO/THTG, PHG, THTM,PHM 

FKCju-FKC’FKC 

FKCT  *  80.645 1/FKCSQ 

SINTHG*SIN(TH0/57.29578) 

COST HG=COS(THO/57. 29578) 

SINPHG=SIN(PHI0/57.29578) 

CO SPHG*C OS (PHI  0/57.295 78 ) 

CALL  GEOTOM(SINTHG,COSTHG,SINPHG,COSPtlG) 
THO-THTM/57. 29578 
PHI0-PHM/57. 29578 

***THO(PHIO  ARE  NOW  GEOMAGNETIC  COORDINATES 
SINPM=SIN(PHIO) 

COSPMsCOS(PHIO) 

SINTN«SIN(TH0» 

COSTH*COS(THO> 

Cl»AU*C0STM»C0SPM*A21*C0STM*SINPN-A3t*SINTM 

C2=A12*C0STM*C0SPM*A22*C0STH*SINPN>A32*SINTM 

C3«A13*COSTM«COSPM»A2J*COSTM*SINPM-AJ3*SINTM 


Ol=-A11*SINPM*A21*COSPM 
D2--A12*SINPM«'A22*C0SpM 
03--  A13*SINPM+A23*C0SPM 

THGTHM=C0STHG*C0$PHG*C1nC0STHG*S1NPHG*C2-SINTHG*C3 

PHGTHH=-SINPHG*C1»C0SPHG*C2 

THGPHM=:C0STHG*COSPHG*01ECOSTHG*SINPHG*O2-SINTHG*D3 

PHGPHM--SINPHG*01+C0SPHG*02 

PR0=SlN(ALPHA/57. 20473) 

PTHG0  =  COS(ALPHA/57.2<!S78)*SIN(BETA/5  7.29578> 

PPGQ  -CO'?  (ALPHA/57.29578)  *COS <HET A/57 . 2957 « ) 

PTHO  -THGTHM*PTHGO+PHGTHM*PPGO 
PP0  =  THG°HM*PTmGG»PHGPHM*PPGO 
RETURN 
END 


SUBROUTINF  INDEX  OF  REFRACTION  ( H , TH , PH I , P P , PT H , PP ) 
COMMON'COORO/THTG.PHG, thtm.phm 
COMMON/INFO/FKCT  , F«C , S IGN 

COMNCM/INFOl/R ,COSTH ,SINTH ,C OSPM ,S INPM , T HO , PH 1 0 , OP T ,S CP  I 

COMMON /INF 02/Y ,YR,YTH, YSQ, NR ,MTH ,X ,XP,0MX,0MX2,XCMX 

COMNON/INF03/MUCK,PRN,PTHN,PPN,COSPS,COSP2,SINP2 ,SINPS 

C0MM0N/INF04/YL, YT.YL2 ,YT2 ,S,C,MUSO,MU 

REAL  NR, NTH, MUCK, HUSO, M'J, N 

R-H*  63  71.0 

COST  H=COS ( TH ) 

SINTH=SIN(TH) 

COSP  M -CO  S ( PH  I ) 

sinpm=sin(phii 

CALL  NAGTOG(SINTH,COSTH,SINPM,COSPN) 

C*  ******  *T  HO ,  PHID  ARE  GEOGRAPHIC  COORDINATES 
TH0=THTG 
PHIQ=PHG 

OPT  = 1 ,  ♦  3.*  C  0STH**2 
SOPT  =5QRT (OPT) 

Y  =  e3'3.*(SOpT*(e37n./R)  ** 3)  /FKC 
C********Y  IS  THE  NORMALIZED  MAGNETIC  FIELD 

YR-  -  3. Q*  Y/R 
YTH=-3.*COSTH*SINTH*Y/CPT 

C********YR,YTH  ARE  DERIVATIVES  WRT  R , TH ,  RESPECTIVELY 
YSO  =  Y*  Y 
NR-  2.*C0STH/S0PT 

C*»*****»Y»NR  IS  THE  MAGNETIC  FIELO  COMPONENT  IN  THE  R  DIRECTION 
HTH-SINTH/SOPT 

C**»*****Y*MTH  IS  THE  MAGNETIC  FIELO  COMPONENT  IN  THE  TH  OT  FRCT  T ON 
C**»*****HE  NOR  CALCULATE  ELECTRON  DENSITY  AND  ITS  DERIVATIVES 
CALL  ELECTRON  DENSITY  <H,N,0NR) 

40  X  =  N*  FKCT 
XR  -0  NR  *F  KC  T 
CMX=1,-X 
30  OMX2  =OMX  *OMX 
XOMX  =X*OHX 

MUCK::SQRT(  PR*PP  +  PTH*PTHNPP*PP) 

C********NOPHALIZE  PR, PTH, PP 
PRNs  PR/MUCK 
PTHNsPTM/MUCX 
PPNrPP/MUCK 
COSPS=PRN*NR»PTHN*HTH 
COSP?=COSPS*COSPS 
SINP2-1. -COSP? 

SI NPS- SORT (SINP2) 

YL  *  Y *C OSPS 

Y  T*Y  *S IMPS 
YL2*Yl»YL 

Y T 2 ‘  YT*YT 
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S=SIGN»SQRT(YT2*YT2«-4.  0*0MX2*YL2) 
0=2.*OMXtS-YT2 
NUSQ=i .0-2.0*XOMX/O 
300  MU^SQRT (MUSQ) 

RETURN 

END 


SUBROUTINE  PARTIAL  MU  (PR,PTH,PP) 

COMM ON/INFOl/R, COSTM,SINTH,COSPM,SINPM,THD, PHI  0, OPT, SOPT 
COMMON/I NFO?/Y,YR, YTH, YSQ,NR,HTH,X,XR,0MX,0MX2 ,XOMX 
COMMON/INF 03 /MUCK, PRN, PTHN ,PPN ,COSPS ,C0SP2 , SIMP2 ,SINPS 
COMM0N/INFO4/YL, YT,YL2,YT2,S,0,MUSQ,MU 
COMMCN/I NF  05/PSPR ,PSPTH, PSPP ,MUPMU  jMURHU ,NUTMU ,MUPHH 
REAL  NR, MTH, MUCK, MUSQ, MU 

REAL  M2YSP.MUXMU, MUYMU ,MUpMU ,HURKU ,MUTHU  »HUPHM 
PSITH=2.*(PR*MTH-PTH*NR)/(MU*SINPS*0PT) 

0SQ=C*0 

M2Y.i'r  =  hUSQ*SINPS 

RMUXM1-2.IH4.0*0MX*YL2/S 

RMUXM2=X0MX*RMUXM1/D 

,*IUX,iU=  (2.0*X-1 ,0-RMUXMZ)  /(  O’MUSQ) 

RMUYMl=2.0*YT2*SINP2+4.0*OMX2*C0SP2 

RMUYM2=-Z. 0*SINP2+RMUYM1/S 

MUYMU=X0MX*Y»RMUYM2/ (OSQ’MUSQ) 

RMUPM1=( YT2-2.0*OMX2)/S-i«0 
MUPMU=XOMX*2.0*YL*YT*RMUPM1/(MUSC*OSQ) 

C  REAPARTIALS  OF  PSI  MRT  PR,PTH,PP 

PSPR-(PR*C0SFS-MU*NR)/M2YSP 
PSPTH=  CPTH*COSPS-MU*MTM) /M2YSP 
PSPP=PP*COSP S/ M2YSP 
C  REASPATIAL  DERIVATIVES  OF  MU 

MURMU=MUXMU*  XK  +rtUYMU*YR 
MUTMU=MUYMU*YTH*MUPMU*PSITH 
MUPHMaO. 0 
RETURN 
END 


SUBROUTINE  ELECTRON  DENSITY  (H,N,ONR> 

COMMON/iCOOP/MGT (1 90 T , EO ( 1 00 ) , KM AX ,HMA X,MTOP 

REAL  N 

HMIN-HGT (1J 

IF(H.GE.HHIN)  GO  TO  1 

N=0.  0 

ONRs  0.0 

RETURN 

C*»**»**»*  INTERPOLATION  BY  PARABOLA 
1  IF(H.GE.HMAY)  GO  TO  4 
IF(H.GE.HTOP)  GO  TO  5 
INTERVAI  SEARCH 
00  10  J3 1 » KM  AX 
XMXI=H-MGT (J) 

XIPHX*HGT ( J* II -H 

IFTTXMXI.GE. 0. 0) .ANO.(XIPHX.GE.0.8)J  GO  TO  3 
10  CONTINUE 
3  I=J 

T1»HGT (I) 

T2*HGT(Iti> 

T33HGT (1*2) 

T12T1J3(T1-T2»*(T1-TJ* 

1 21T  23*  t  T2-T li •( T  2-T3T 
T31T3?M’3-TU  •  ( T3*T 2) 
c. 8m  VECTORS  FOLLOW  P1(P2<P3 
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Pi  =  ( H- T2 ) * (H-T  3  > /T12T13 
P2=(H-T1)*(H-T3)/T21T23 
P3= ( H-Tl ) * (H-T2) /T31T32 

C*********  ELECTRON  DENSITY 

N=£0(I)*P1+EO(I+1)*P2EEO (I  *2 ) •  P3 
N=10Q0.0*N 

c, ELECTRON  DENSITY  DERIVATIVE 
A1=(2.*H-T2-T3)/T12T13 
A2=(2.*H-T1-T3)/T21T23 
A3=(2.»H-T1-T2'/T31T32 
ONR=Al*EO(I) ♦A2*EDCIfl)»A3*En(I»2) 
ONR=  1000  »0*DNR 
RETURN 

5  T1  =  HGT ( I  — 1 ) 

T2SHGTCI) 

T3=HGT (1*1) 

T12T13=(Tl-T2)*(Tl-T.3) 

T21T23=(T2-T1)*(T2-T3) 

T31T32=(T3-Ti)*(T3-T2) 

Pl-(H-T’)* (H-T  3) /T12T13 
P2=(H-T1 )• (H-T  3) /T21T23 
P3=( H- Tl)* (H-T 2) /T 31 T3 2 


c, electron  oensity 

N=E0 (I-l)*Pl+ED( I) »P2*E0f 1+1) *P3 
N=1000  .'J*N 

C****»»**»  ELECTRON  DENSITY  DERIVATIVE 


A1=C2.*H-T2-T3)/T12T13 

A2-<  2. *H -T 1-T3J/T21T23 

A  3  = I  2 . *H-T 1-T2J/T31T32 

ONR=Al*E 0(1-1)  ♦A2*tO(I>»A3*EO(I+l> 

DNR=1000.0»DNR 

RETURN 

C»»*»***»»  H  GT  HN AX 

4  N=ED (KMAX) 'lOOO. 0»CNR* <h -HMSX) 
RETURN 
ENO 


SUBROUTINE  OATA  READ  IN 

CDMMCN/SrOOP/HGT  (100)  ,E(1  (100)  , KM AX ,HH AX, HTOP 

00  10  K- 1  *  10  0 

REAO  (5,1)  HGT(K)  , ED ( K) 

1  FORMAT <FB. 3, 1GX, F10. 4) 

C*» LAST  DATA  CARD  HAS  H=2000 
TF  IH'jT  (Y)  i  CT  ,1  £3  ji  u)  Go  to  lx 

10  CONTINUE 

C*»***»»»«  «MAX  IS  THE  TOTAL  NUMBER  OF  PAIRED  OATA  POINTS 

11  KMAX--K-1 

HMAX  *MGT (KMAX) 

MTOPsHGT (KMAX-1) 

C***»»**»»  VERIFY  OATA  OY  PRINT  OUT 
WRIT  E  ( 6,  2) 

2  FORMAT (1M1 ) 

WRITE ( 6, 3) 

3  FORMAT (2H  , 4X,>HH(KM) ,10x,i2hEL ECTRONS/CC // ) 

DC  20  i , KMAX 

WRITE (6,4)  HGT  (X) ,EO(K) 

20  CONTINUE 

4  FORMAT (F  11. 3,1  OX, 613. 5) 

RETURN 

ENO 


SUBROUTINE  GEOTOH(l(>.TNG,COSTHG,SZ nPHG,COSPMGI 


COMMON/COORD/THTG.PHG, THTM.PHM 

A=COSPHG*SINThG 

8-SI  NPHG’S  IN  i  , -r 

COSTM-0.  07  14822*  - (1 «  18£l278»n<  0.  97992*C0STHG 
SINT  M=S0RT (1 .0-CCSTM*COSTM) 

IP (COSTN.LT. 0 . 0)  GO  TO  5 
THTM=ATAN(SINTn/COSTM) *57. 29  578 
GO  TO  4 

5  THTM-ATAN(SINTM/COSTNI  *8  7.  28578*  110.  0 
4  S.'NPMr  (0. 83358 *A  *0.858  77»9)  /SI NT7 

COSPM=(0. 3511739*8-9.914 «3 3 7*0-0 .1 9937*COSTHG> /SINTM 
IF (GOSPN.LT. 0. 0>  GO  TO  6 
PHM=  4T  AN (SINPH/COSPN)* 57 .795  39 
GO  TO  7 

6  PHMrATAM  (SINRM/COSPm>*57.?9579*180.0 

7  IF (°HM.L  T. 0 . 0)  PHM^PHM ♦ 3  60 . U 
RETURN 

ENO 

EOF 


SUBROUTINE  SFTPOT  (A;)OR,COEF) 
EXTENDEO  AOOR 

IF(COrF.FO.l.uOOQ)  COE  F  =  0 , 99  99 
IF  (OCEF.CQ-O.OOUO)  OOCF-rQ.  0002 
ITRY  =9 

2  ITPY=ITRY*1 

IF(ITRY.GT.3)GOT04 

CALL  OSTPTO<1,AOCR,COFF,IERR) 

IF(TERR.FO.l)  PETURN 

IF ( IERR. E0.2 )  NR  ITE  (6,  20 )  AOOR 

IF  (IERR.EQ.2)  PAUSE 

IF  (IfRR.EQ.3)  WRITE (6, 30)  AOOR 

IF ( IERR. EQ . 3 )  PAUSE 

IF  (IERR.EQ.4)  WRITF (6  >40)  AOOR 

If  (IERR.EQ.4)  PAUSE 

IF (I£RR.EG.5)GOTC2 

RETURN 

4  WRITE ( 6( 50 ) AOOR 

IF  (IFRR.EQ.5)  PAUSE 
20  FORMAT  ( 1 7H  INVAlIC  A'jOPCSS  ,A4) 

30  cORM  A  Y ( 1 6H  CCEFF  CVtRFLOW  ,A4) 

40  FORMAT  ( 21 H  CONSOLE  01  SO ONNECT EO > 
50  FORMAT  U4H  NULL  FAILURE  ,A 4) 
RETURN 
END 


SUBROUTINE  MAGTOCC^INTHM^OSTHM,  SI NPHM .COSPMM' 
COMMON/O 0090/ T  HTG.PMG, IhIM ,PHM 
A*SINTHM*COSPMF 
0sSINTM9*SINFHN 

COST G- 9. 97 99 2* COST  MM -0.19937* A 
SINTG^SlRf  <1 .0 -CCSTG*OCSTG) 

IF (COSTS. LT. 0. 0)  GO  TO  15 
ThTG =A TAN (SI  NT  0/ COSTS) *5  7.295/8 
CO  TO  14 

J  5  ThTG * A  TAN (SI  NIG/ COST G) *5  7.29578*  190.  0 
14  S I NPC®  (-0.9148337*8*0.  .(58T7*fl"0.  186  J  2  78*  COS  T  MM ) /S I NTG 
COSPG*  (0 . 351  17  39*A*0.91l58*iT*a.0  7149  22*COSTMM)/SlNTG 
IF(CC3PC..LT.O.O>  GO  TO  IE 
P*)G=»TAn  (S INPG /C  CSPO ) *  57 . 7  95  7  i 
GO  TO  17 

If  Pm5*ATAN(SInPG/CCSP.)»57.?957i«J 80.3 
17  IF (PMG  .L T. 0. 0)  Pw6 iPNG  *360.8 
Rt  TURN 
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APPENDIX  II  B 


IIB1.  PROCEDURE  FOR  RUNNING  THE  RAY  TRACING  PROGRAM  (I) 

i)  Load  main  program  plus  subroutines  and  data  into  the  hopper. 

ii)  Type  $G0  on  the  console  typewriter. 

Sequence  of  events 

a)  Program  plus  data  will  be  read  in. 

b)  After  reading  data,  the  electron  density  profile 
used  will  be  printed  plus  the  target  and  trans¬ 
mitter  parameters. 

c)  Next  the  run  number  will,  be  printed  with  the 
initial  conditions  which  specify  an  ordinary 
or  an  extraordinary  ray. 

d)  Tlis  should  be  followed  by  the  IC  mode  of  the 
analog  computer  and  the  console  typewriter 
output  message:  "BLIP  FSW  1012  FOR  IC  PRINT... 

PRESS  FLAG  8".  This  will  be  followed  by  a 
PAUSE . 

iii)  Release  the  PAUSE.  There  are  two  courses  of  action  which 

could  follow: 

a)  If  FSW  1012  was  "blipped"  an  IC  print  will  follow 
ending  with  a  PAUSE.  Releasing  this  PAUSE  will 
cause  the  typewriter  to  type,  "PRESS  FLAG  8  TC 
CONTINUE". 

b)  If  FSW  1012  was  not  "blipped"  the  typewriter  will 
immediately  type,  "PRESS  FLAG  8  TO  CONTINUE" 
after  releasing  the  PAUSE  of  step  ii)J. 

c)  Either  (a)  or  (b)  of  this  section  is  followed  by 
a  PAUSE. 
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iv)  Again  release  the  PAUSE. 

Sequence  of  events 

a)  Once  the  toleranqe  on  PHI  is  satisfied,  i.e.  the 
end  of  one  run,  the  analog  computer  will  be  placed 
into  HOLD. 

b)  The  line  printer  will  write  the  headers  for 

"I,  ALPHA,  BETA,  H,  THD,  PHID"  followed  by  their 
values  at  the  end  of  the  run. 

c)  Next  the  IC  potentiometers  will  be  set  (H,  PHID, 

THD),  the  analog  computer  will  go  to  IC  and  the 
typewriter  message  of  step  ii)d  will  be  repeated. 

Also,  as  a  verification  of  the  next  £,  Ot  pair 
their  modified  values  will  be  printed  on  the  line 
printer. 

v)  From  this  point  steps  ii)d  through  iv)  are  repeated  until 
convergence  criteria  are  met. 

vi)  When  the  ray  is  the  optimum  one  within  the  constraints,  the 
typewriter  output  is  "THIS  IS  THE  OPTIMUM  RAY  PATH".  A 
PAUSE  will  follow. 

vii)  Release  the  PAUSE  and  repeat  steps  ii)d  through  Iv)  for 
this  ray.  (We  are  repeating  the  optimum  ray  for  time 
history  storage.) 

a)  At  the  end  of  this  ray  path  the  typewriter  message 
is,  "SET  SLO  FOR  TIME  HOSTORY  PRINT"  followed  by  a 
PAUSE.* 

viii)  Set  FSW  1012  and  release  the  PAUSE, 

a)  Time  history  printout  will  follow, 

b)  Typewriter  message  at  conclusion  of  print  is, 

"TIME  HISTORY  PRINTED  RESET  SLO".  (The  last 
message  is  to  ensure  that  th<  switch  is  not 
left  in  the  set  position  for  the  next  run.) 


ix)  This  is  conclusion  of  one  run  for  a  given  electron  density 
profile, 

x)  To  load  new  electron  density  data,  replace  the  old 

electron  density  profile  with  a  new  one  in  the  data  deck. 
The  end  of  electron  denfdty  data  signal  (H  *  2000)  must 
be  retained. 

*N0TE:  Since  line  zero  (SLO)  is  FSW  1012. 


IIB2. 


CHANGES  TO  ANALOG  BOARD  REQUIRED 
TO  CONVERT  FROM  AUTOMATIC 
TO  MANUAL  RAY  TRACING 


Amp  210  goes  directly  into  trunk  line  330  for  the  manual 
case.  In  the  automatic  case  amp  414  goes  into  trunk  330  and  amp 
210’s  output  into  trunk  330  in  removed. 

For  the  automatic  case,  the  output  of  amp  014  should  input 
comparator  000. 

For  the  manual  case,  the  output  of  amp  812  should  input 
comparator  000. 

For  both  cases  the  timer  in  the  8400  should  be  set  at  10  usee. 

Other  than  the  above  there  is  no  change  to  the  analog  board.  The 
static  test  program  is  applicable  to  both  cases.  Essentially  all 
that  is  done  by  the  above  is  to  remove  the  H  optimum  circuit. 

The  pots  will  be  set  by  the  static  test  program  but  it  doesn't  make 
any  difference  since  amp  210  bypasses  the  optimum  circuit. 


IIB3. 


MANUAL  RAY  TRACE 
PHASE  II 

SWITCH  AND  SENSE  LINE 
ASSIGNMENT 


Sense  Line  F.  Switch 


Function 


0 

1012 

*1 

1211(1) 

*2 

1053(0) 

3 

1013 

4 

Not  Used 

5 

l  ot  Used 

*6 

OP 

*7 

IC 

IC  Printout  Check 
Pj  *  f (Pr j,  PQ,  n) 
P0  *  f(Pr,  P$,  u) 
New  Run 


Hold  Simulated 
IC  Simulated 
L-Alt.  C -Phase  Path 
R-Group  Path 


411 

(For  X-Y  Plotter) 

♦Asterisks  denote  hard  wired  senee  lines 
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IIB4. 


PROCEDURE  II 

(ORIGINAL  MANUAL  RAY  TRACING) 


i)  Place  deck  (Main  plus  subroutines)  in  the  hopper 

ii)  $G0  on  console  typewriter 
Sequence  of  Events 

a)  Program  plus  data  will  be  read  in 

b)  After  reading  data  (see  data  deck  org) 

1)  The  electron  density  profile  will  be  output 
on  the  line  printer 

2)  The  IC  pots  (C310,  C701,  C910)  will  be  set 

3)  The  analog  computer  will  go  to  the  I.C.  mode 

iii)  To  make  sure  the  program  has  been  loaded  properly,  blip 
function  switch  1012.  This  will  give  you  an  IC  printout 
followed  by  a  Fortran  pause  (Flag  8  high).  Release  the 
pause,  more  printout  will  insue  with  a  second  pause. 

Release  the  second  pause. 

iv)  Go  to  the  analog  console  and  manually  place  the  analog 
computer  in  the  operate  mode.  The  ray  will  run  and  the 
analog  computer  will  hold  at  the  distance  607  KM.  i.e, 

is  distance  from  Green belt.  MD.  to  AFCRL,  Bedford,  Mass. 

If  at  this  point  a  printout  is  desired  follow  the  stens  in 
iii  After  this  place  the  analog  computer  into  IC.  This 
is  the  end  of  on®  ray. 

v)  To  reinitialize  the  program  for  another  ray  (i.e.  increase 

alpha)  Flip  function  switch  1013,  to  the  left.  The  new  data 
cards  will  be  read,  their  values  printed,  the  pots  set  and  'che 
analog  will  go  to  IC.  From  here  repeat  step  iv) . 


II B5. 


DATA  DECK  ORGANIZATION 


Data  is  read  into  the  computer  at  the  beginning  of  each  run. 
Initial  conditions  for  each  run  are  contained  on  three  cards 
as  indicated  below: 

CARD  (1):  .RUN  (18) 

I  TUN  is  the  run  number.  It  is  fixed  point  and 
can  be  up  to  eight  digits  in  length. 

CARD  (2):  h,  THO ,  PHIO,  ALPHA,  BETA  (5F1S.7) 

The  initial  value  of  Height,  Theta,  Phi, 

Alpha  and  Beta  are  on  card  2  in  floating 
point  form.  They  can  be  defined  with  up 
to  seven  decimal  places. 

CARD  (3):  SIGN,  FKC  (2F15.7) 

SIGN  indicates  the  type  of  ray.  For  an 
ORDINARY  ray,  SIGN  =1.0 

For  an  EXTRAORDINARY  ray,  SIGN  =-1.0. 

FKC  is  the  frequency  of  the  ray  in  kilohertz. 
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I II. 1.0  INTRODUCTION 


A  nuclear  emulsion  is  a  material  which  records,  photographically,  the  tracks 
of  charged  particles.  An  ionizing  particle,  one  with  sufficient  energy,  on 
encountering  a  crystal  of  nuclear  emulsion  renders  if  developable.  After 
development  and  further  processing,  the  paths  of  charged  particles  that 
penetrated  the  emulsion  are  visible  through  a  microscope  as  trails  of  minute 
grains.  This  trail  of  grains  represents  a  three-dimensional  image  of  the 
charged  particle's  path. 

The  major  application  of  nuclear  emulsion  is  in  experimental  physics. 
Emulsions  provide  the  means  with  vshich  interactions  between  charged  parti¬ 
cles  may  be  observed. 

The  instrument  used  for  analyzing  the  behavior  of  tracks  in  emulsions  is 
the  microscope.  The  path  of  a  high  energy  particle  through  an  emulsion  is 
presently  scanned  by  human  scanners  using  microscopes.  This  scanning 
process,  as  performed  by  humans  is  tedious  and  subject  to  errors  caused 
by  fatigue  of  the  human  scanner. 

This  report  describes  the  concepts  developed  for  a  track  tracing  system 
and  their  embodiement  within  an  optically  scaled  breadboard  model  for 
automating  the  scanning  and  analysis  of  nuclear  emulsions.  The  results 
obtained  with  this  breadboard  model  established  the  credence  of  these 
concepts . 
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III. 2.0  REQUIREMENTS  OF  AN  AUTOMATED  SCANNING  SYSTEM 

Any  system  for  automating  the  scanning  process  for  nuclear  emulsions 
must  coordinate  and  integrate  the  following  tasks: 

a)  Search  the  emulsion  for  track  entries; 

b)  Determine  whether  the  entry  point  is  isolated  or  part 
of  a  track; 

c)  Determine  directional  properties  of  tracks  emanating 
from  an  entry  point; 

d)  Trace  along  tracks  until  a  termination  or  vertex  is 
encountered.  A  vertex  is  defined  as  a  point  along  a 
track  where  splitting  occurs;  this  also  includes 
degenerate  splitting,  i.e.  no  splits  at  all. 

e)  Determine  whether  a  termination  along  a  track  is  either 
an  end  point  of  the  track  or  a  vertex;  if  a  vertex  is 
detected,  the  emanating  tracks  must  be  classified  accord¬ 
ing  to  their  directional  properties  for  subsequent  track 
tracing; 

f)  Organize  detected  entry  points,  tracks,  vertices,  and 
terminations  in  a  storage  format  which  enables  the 
assemblage  of  several  frames.  Nuclear  emulsions  are 
generally  stacked  in  frames,  frames  are  examined  indi¬ 
vidually,  so  that  the  information  retrieved  from  a  frame 
must  be  coordinated  with  information  retrieved  from 


upper  and  lower  frames. 


Figure  (I1I-1)  depicts  a  conceptual  block  diagram  of  the  optically  scaled 
breadboard  model . 


FRAME 


FIGURE  (III-J).  CONCEPTUALIZATION  OF  TRACK  TRACING 

The  optical  system  comprises  the  collection  of  lensei  which  project  real 
images  of  events,  from  the  nuclear  emulsion  onto  the  face  of  an  optical 
sensor.  These  images  must  represent  with  minimal  distortion,  the  real 
events  within  the  emulsion.  The  breadboard  model  used  a  simple  Cooke 
triplet  with  1:1  magnification  and  a  .375  inch  depth  of  focus. 

The  optical  sensor  produces  electronic  signals  from  images,  projected 
onto  its  face,  of  events  or  segments  of  events  within  the  nuclear  emulsion. 
These  electronic  signals  represent  the  observed  measurements  taken  from 
the  nuclear  emulsion  --  figuratively,  it  is  the  "eye"  of  the  track  tracing 
system.  Our  breadboard  simulation  used  an  electiostatically  focused 
viJicon  as  the  optical  sensor. 
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The  means  for  controlling  the  relative  movement  of  the  frame  with  respect 
to  the  optical  system  and  sensor  as  well  as  for  movement  along  an  optical 
axis  (i.e.  an  axis  along  which  optical  properties  are  varied  for  depth 
perception)  was  accomplished  by  servomechanisms. 

A  hybrid  computer  was  selected  as  the  medium  for  implementing  the  bread¬ 
board  model  of  the  system  because  it  combined  continuous  and  sequential 
operations  along  with  storage  of  data  --  precisely  the  requirements  of 
an  emulsion  scanning  system. 

The  analog  sections  provide  for  the  generation  of  scanning  regimes,  data 
retrieval  from  the  vidicon,  and  generation  of  control  functions  for  the 
servomechani sms . 

The  digital  section  serves  as  an  executive  to  coordinate  the  scanning 
tasks,  to  provide  for  information  retrieval,  and  its  transformation  to 
a  data  base  compatible  with  the  assemblage  of  data  from  the  many  frames 
of  a  stack. 

The  implementation  for  the  tasks  listed  at  the  beginning  of  this  section, 
were  combined  within  our  optically-scaled  breadboard  model  for  a  nuclear 
emulsion  track  tracing  s_>  -tern.  The  remainder  of  this  report  treats  the 
scanning  concepts  applied  to  this  model  and  their  implementation  and  the 
implementation  of  a  data  base  for  the  assemblage  of  data  from  frames  of 
a  stack. 
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Ill, 3.0  CHARACTERISTICS  OF  THE  OPTICALLY  SCALED  BREADBOARD  MODEL 


During  the  early  phase  of  this  study  numerous  track  following  and  vertex 
analysis  procedures  were  studied  by  EAI  in  an  effort  to  select  methods 
suitable  for  the  analysis  of  nuclear  emulsion  stacks  by  hybrid  computer  -- 
vidicon  sensor  systems.  The  procedures  ultimately  selected  and  incor¬ 
porated  in  the  logic  of  software  programs  subsequently  developed  are 
based  upon  a  "vertex -to-vertex1'  philosophy,  administrative  control  of 
which  resides  in  a  FORTRAN  IV  main  program.  Three  basic  scanning  modes 
are  incorporated .  They  are: 

Edge  Scanning 

For  purposes  of  this  study,  an  edge  was  defined  as  the  physical 
boundary  of  the  frame.  Track  intersections  with  an  edge  are 
located  by  optically  scanning  a  rect;  gul&r  grid  imposed  over 
the  edge  surface. 

Vertex  Analysis 

Vertices  thought  t>  ist  at  a  point  '’re  analyzed  by  optically 
constructing  a  thin  spherical  she  1  about  the  point  in  question, 
Confirmation  of  cracks  suggested  by  "blobs"  encountered  within 
the  spherical  shell  itself  is  accomplished  by  scanning  along  a 
line  connecting  such  blobs  with  the  central  point. 


Tracx  Following 


A  track  is  followed  from  a  starting  point  with  an  initial 
direction  by  scanning  along  three  directions  --  the  initis* 


direction  and  two  adjacent  ones  --  selecting  as  a  valid  direc¬ 
tion  that  angle  for  which  a  maximum  optical  length  occurs. 

This  length  and  angle  determine  the  next  position  along  the 
track. 

The  administrative  control  of  the  scanning  procedures,  recording  of  data 
and  its  correlation  from  a  stack,  consisting  of  an  arbitrary  number  of 
frames,  is  accomplished  by  the  main  program. 

III. 4.0  THE  "VERTEX -TO-VERTEX"  PHILOSOPHY 

An  event,  as  observed  in  nuclear  emulsions,  is  a  collection  of  verticles 
and  links.  These  links,  in  general,  are  curved  lines  as  opposed  to 
straight  lines,  so  that  a  characterization  of  an  event  should  include 
vertices  and  links,  coupled  with  a  measure  of  curvature  for  the  links 
connecting  the  vertices. 

The  events,  based  on  the  characterization  above,  can  be  recorded  in  a 
compact  form  which  utilizes  the  vertex  locations  and  their  multiplicity 
along  with  the  directional  properties  of  the  links  connecting  the 
vertices . 

This  characterization  then  defines  sequences  within  the  track  tracing 
process,  namely: 

1)  Detection  oi'  track  entries  --  treated  as  vertices. 

2}  The  directional  properties  of  links  emanating  from  a  vertex.  The 
number  of  emanating  links  is  termed  the  multiplicity  of  the  vertex. 


3)  The  tracing  of  links  emanating  from  one  vertex  to  another  vertex, 
which  may  have  multiplicity  zero  (a  termination  point,  or  a 
multiplicity  of  two  or  greater  (a  splitting  of  the  track  indicating 
disintegration  phenomena). 

These  regimes  were  incorporated  into  the  software  for  the  optically  scaled 
breadboard  model.  The  administrati  e  control  for  coordinating  the  inter¬ 
action  of  these  regimes  resides  within  the  main  program.  These  regimes 
collectively  comprise  what  we  term  the  "vertex -to-vertex"  philosophy  of 
the  breadboard  model. 

1 1 1. 5.0  THE  VIDICON-OPTICS  SU3SYSTEM 

The  vidicon  is  a  pbotoconductive  device  which  generates  electrical  current 
proportional  tc  the  light  intensity  incident  or.  its  face  --  "  photoconduc- 
tive  layer  of  phosphor.  An  image  projected  into  the  face  of  the  vidicon 
can  be  transmitted  electrically  by  exploring  the  face  in  a  systematic 
manner  and  transmitting  at  each  instant  the  generated  current.  The  result 
of  such  a  process  is  to  produce  a  current  that  varies  with  time  in  accord¬ 
ance  with  the  light  intensity  of  successive  elements  of  the  image  projected 
into  the  vidicon's  face.  This  process  of  exploring  an  image  to  obtain  a 
current  that  varies  with  time  in  accordance  with  the  light  intensity  of 
successive  areas  of  the  image  is  called  scanning. 

The  scanning  „f  the  vidicon's  face  is  accomplished  by  sweeping  an  electron 
beam  across  the  nhotoconductive  layer  of  the  vidicon's  face  by  means  of 
electrostatic  deflection  voltages. 


Schematically,  the  vidicon  face  is  represented  below 


VOLTAGE 


The  deflection  voltages  control  the  position  of  the  electron  beamts  point 
of  impact  with  the  photoconducti ve  layer.  The  current  produced  by  the 
photoconductive  layer  is  proportional  to  the  light  intensity  incident  in 
the  layer  at  the  point  of  impact  for  the  electron  beam. 


The  scanning  mode  used  for  the  breadboard  model  vidicon  ’ s  a  radial  scan 


with  a  variable  radius  and  variable  angle.  Schematically, 


The  angle  a  is  measured  relative  to  a  horizontal  axis.  This  scanning 
procedure  was  realized  by  generating  the  x  and  y  deflection  voltages 
in  the  form: 

^volts  Ry)  COSa 

^volts  =  *  Rv)  s^na 

with  Rv  voltage  represe^-ing  the  radius  of  the  vidicon  face,  and  R  a 
swept  voltage  varying  from  zero  to  2Ry  every  2  milliseconds. 

The  analog  implementation  of  this  scanning  procedure  appears  in 
Figure  (1 1 1-2) . 


FIGURE  (1 1 1 -2).  RADIAL  SCAN  VOLTAGE  GENERATION 


The  waveform  of  R  produced  by  the  above  oscillator  is  shown  in  Figure  (III-3) 


FIGURE  (1 1 1 -3).  WAVE  3HADE  G£  THE  RADIAL  SCAN  VOLTAGE 

Information  is  retrieved  during  the  rising  slope  of  the  radial  voltage, 
on  the  descending  slope  the  vidicon  is  blanked  to  reduce  the  average 
current  produced  by  the  vidicon  --  a  protective  measure. 

The  Vidicon  Output 


We  consider  the  vidicon  face  to  be  a  disk  onto  which  images  of  segments 


the  disk  fixes  a  reference  frame  for  projected  images  of  events. 


For  each  scan  angle,  a,  the  vidicon  scan  of  its  face  --  for  the  projected 
image  above  --  produces  a  current  proportional  to  the  intensity  of  light 
incident  on  its  face.  Typical  vidicon  outputs  are  shown  below: 


(l)o*0  Y 


VIDICON  SCAN  ENCOUNTERS  A 

CURRENT  /LINE  OF  THE  EVENT 
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REGION  REGION 

X  IS  DIRECTLY  PROPORTIONAL  TO  TIME 


VIDICON  SCAN  ENCOUNTERS  A 

CURRENT  /LINE  OF  THE  EVENT 


X*— R  v  Y=Y*  X=Rv 


Y  IS  DIRECTLY  PROPORTIONAL  TO  TIME 
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Examination  of  the  above  scan  shows  that  when  the  scan  encounters  a  dark 
region,  such  as  a  portion  of  a  line,  the  light  intensity  decreases  causing 
the  current  to  drop  to  zero  producing  a  pulse  shaped  output.  The  location 
of  these  pulses  relative  to  the  origin  of  the  established  coordinates  is 
determined  from  the  time  at  which  the  scan  encountered  t lie  time  segment. 

I II. 6.0  EDGE  SCANNING  IN  THE  BREADBOARD  MODEL 

Our  study  assumed  that  the  stack  of  frames  was  shielded  from  above  and 
below  and  that  all  observable  events  begin  on  the  sides  of  the  stack. 

This  assumption  covers  the  largest  class  of  expected  events.  These  are, 
however,  events  which  are  visible  only  within  the  interior  of  a  stack, 
o.g.  non-interacting  primary  particles  which  decay  within  the  stack  and 
produce  interacting  particles,  that  is,  non-interacting  particles  are 
not  visible  as  tracks. 

With  this  assumption  that  events  begin  on  the  edge  of  a  frame  in  the 
stack  the-  track  entry  point  is  located  by  imposing  a  rectangular  grid 
of  cells  over  the  edge  of  a  frame  ar.d  determining  which  it'lls  have  been 
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erposed.  That  is,  the  entry  point  is  not  really  a  point  but  a  region  of 
exposed  grams  of  emulsion  --  a  blob  --  on  the  side  of  a  frame.  The 
exposed  region  must  be  enclosed  by  a  closed  curve  and  its  centroid  used 
as  the  point  of  entry  of  a  track. 


FIGURE  (1 1 1 -4 ) .  THE  BLOB  CONFIGURATION  ON  THE  EDGE  OF  A  FRAME 

Figure  (IIJ-4)  depicts  the  edge  of  a  frame  with  a  rectangular  grid  imposed. 
The  rectangular  grid  is  generated  by  the  points: 


-  iAx 

i  =  0,1,2... 

■  »N 

=  iAz 

i  =  0,1,2,.. 

.,M 

The  z^  noints  refer  to  the  plane  of  focus  of  the  optical  system  with  Az  its 
depth  of  focus .  Ay  correspcnds  to  the  width  of  the  electron  beam.  Ax  is 
the  scan  window  for  the  vidicon  face,  e.g.,  Ax  =  2Ry  for  a  coarse  grid, 

Ax  =  K2Ry,  K<1  for  a  finer  grid. 


120 


A  blob  is  indicated  by  the  darkened  cells  within  the  grid.  Generally, 
the  blobs  do  not  occupy  cells  in  so  orderly  a  manner,  but  for  our 
purposes  this  sketch  does  not  impose  any  difficulties. 

For  the  breadboard  model,  the  track  entries  were  represented  as  the 
intersection  of  line  segments  from  events  and  the  physical  boundary  of 
the  frame  face. 


F 1  CURE  (111-5).  IMF.  BR  HAL)  BOARD  MOD  ML  IMPLEMENTATION  OF  A  FRAME 

Figure  (111-5)  depicts  the  implementation  of  a  frame  in  the  breadboard  model. 
Three  levels  of  z  were  u-^ed  --  due  to  limitations  in  physical  site  and 
optical  parameters  of  our  components.  Each  level,  or  2 -phase  of  a  frame 
comprised  a  view  graph  slide  with  events  represented  by  thin  strips  of 
tape.  Three  dimensional  tracks  were  represented  by  taping  projections 
in  each  2- level . 


The  edge  scanning  proceeded  by  moving  the  frame  relative  to  the  vidicon 
optics  assembly. 


FIGURE  (II I -6).  THE  EDGE  SCAN  TECHNIQUE 

The  rectangular  shape  in  Figure  (I I I -6)  represents  the  physical  boundaries 
of  a  frame.  The  circles  imposed  along  these  edges  represent  positions  of 
the  vidicon  optics  assembly  relative  to  the  frames  edges.  Ihe  edges  of 
the  frame  are  numbered  1,2,3,  and  4;  respectively;  for 
x  a  0,  0ixI*max;  x  *  xmax >  0<y<ymax;  y  »  ymaLX,  0<x<xmax;  x  =  0,  0<y<ymax. 

This  edge  scan  program  is  illustrated  by  Figure  (I II -7) 
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FIGURE  (1 1 1-7) .  ILLUSTRATION  OF  EDGE  SCAN  PRODEDURE 


The  electron  beam  when  scanned  across  the  vidicon  face  with  a=0  produces 
a  pulse  on  encountering  the  line  segment. 

pulse  indicating  intersection  of  the 
scanning  beam  with  the  projected  image 
of  the  line  segment. 

The  time  of  occurrence  of  this  pulse,  relative  to  the  start,  determines 
the  coordinate  of  this  entry  point  in  the  reference  frame  of  the  vidicon 
face.  Upon  detection  of  an  entry  point,  the  main  program  proceeds  to  a 
vertex  analysis.  The  vertex  analysis  determines  whetner  or  not  tracks 
emanate  from  the  point  and  if  so,  their  number  and  respective  direction 
angles . 

1 1 1.7.0  VLRTEX  \N A LYSIS  IN  THE  BREADBOARD  MODEL 

Once  a  vertex,  such  as  a  track  entry  point,  is  detected,  the  vertex  is 
analyzed,  or  processed,  to  identify  emanating  lines. 

The  process  begins  by  positioning  the  center  ot  the  vidicon  optics  assembly 
over  the  vertex  point. 


Assuming  an  origin  0(x=0,  y=0,  z=0),  the  first  cell  scanned  is  cell  1 
(zj,  0<_x<xi),  if  a  track  entry  is  detected,  we  proceed  to  a  vertex 
analysis  of  the  entry  point,  if  not,  we  scan  cell  2 (z-? ,0<x<Xj)  by  moving 
the  vidicon  optics  assembly  along  the  z  axis,  if  no  track  entries  are 
detected  we  scan  cell  3 (zj.Ofx^xj) ,  again  if  no  track  entry  appears  we 
scan  cell  4(Z3,xj<x<X2)  by  moving  the  vidicon  assembly  to  a  n ew  starting 
position  along  the  frame  edge  (x=x^),  and  proceed  to  scan,  provided  no 
entries  are  detected,  cells  5,6,...,  etc. 

Detection  of  an  Entry  Point 

An  entry  punt  is  defined  the  intersection  of  the  frame  boundary  and 
a  line  segment  of  an  evert.  As  we  scan  along  the  edge  of  a  frame,  as 
shown  above,  images  are  projected  onto  the  vidicon  face. 

Y I 

LINE  SEGMENT 
FROM  EVENT 

- - X 

E08E  I 

The  first  position, of  the  assembly,  circle  1,  will  not  have  a  line  segment 
projected  onto  its  face.  The  second  position,  circle  2,  of  the  assembly 
will  have  the  image  of  the  line  segment  projected  onto  its  face. 

a  *0 


An  image  of  the  vertex  point  and  its  neighboring  area  is  projected  onto  the 
face  of  the  vidicon  optics  assembly.  The  vidicon  face  is  scanned  for  suc¬ 
cessive  a  scan  angles,  for  each  such  scan  angles,  for  each  such  scan  angle, 
an  integral  derived  from  the  vidicon  signal  is  computer,  referred  to  as  the 
Brightness  Integral. 


Vidicon 

Signal 
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With  a  vidicon  signal  as  sheen  above,  a  logic  signal  is  derived  by  com¬ 
parison  against  a  threshold,  Turing  the  data  retrieval  portion  of  the 
radial  scan.  This  logic  signal  controls  the  operate  mode  of  an  integrator 
with  constant  input.  So  that  the  integrator  value,  at  the  end  of  a  scan, 
represents  a  measure  of  coincidence  for  the  scan  with  a  line  segment 
projected  on  the  vidicon  face. 

I f  we  graphically  represent  the  results  of  this  procedure,  the  brightness 
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Peaks  occur  at  the  scan  angles  a ^  and  along  which  the  projected  image 
has  line  segments.  The  brightness  integnl  also  serves  to  establish, 
when  compared  against  a  threshold  value,  the  existence  of  a  line  segment. 
These  peaks  are  detected  by  the  main  program  --  at  present  the  program 
.  s  capable  of  recognizing  as  many  as  five  peaks  --  but  this  can  be  extended 
to  include  greater  numbers. 

Each  of  the  detected  line  segments  emanating  from  a  vertex  print  is  then 
traced  along  its  path  until  a  new  vertex  —  either  a  termination,  or  a 
splitting  point  --  appears,  at  which  time  the  vertex  analysis  procedure 
is  repeated. 

III. 8.0  TRACK  TRACING  ON  THE  BREADBOARD  MODEL 

The  track  tracing  regime  occurs  after  the  disclosure  of  non-zero  direction 
angles  for  emanating  tracks  of  a  vertex.  The  track  trace  is  initialized 
with  a  vertex  point  and  a  direction  angle.  The  track  tracir-"  continues 
until  a  new  vertex  is  encountered. 

In  the  diagram  below,  A  is  the  vertex  point  --  the  starting  point  of  the 
track  tracing  process.  0  is  the  center  of  the  vidicon. 


'Hie  vidicon  optics  assembly  is  positioned  at  the  point  0,  computing  +he 
assembly's  coordinates  as 


xv  =  XA  +  Rv  cosct 

yv  =  yA  +  Rv  sinct 

The  heavy  dashed  line  in  the  diagram  represents  the  projected  ir  age  of  a 
line  segment  emanating  from  the  vertex  of  point  A  with  direction  angle  a 

The  vidicon  scans  from  point  A  along  the  angle  u.  The  corresponding 
vidicon  signal  appears  as  a  series  of  pulses,  each  pulse  representing  an 
encounter  of  the  electron  beam  with  a  dashed  section  of  the  line. 


Vidicon 


A  0 

- ►  SCAN  LENGTH 

That  scan  length  for  whcih  pulses  no  longer  appear  is  recorded  as  the 
maximum  travel,  Rscan,  from  point  A  along  the  direction  angle  a  of 
coincidence  between  the  line  segment  and  scan.  For  the  case,  shown  in 
the  above  figure,  Rscan  =  2Ry. 

For  each  starring  point  of  a  known  line  segment,  such  as  A,  three  suc¬ 
cessive  angular  scans  take  place;  w-Aa,  a,  a+Aa  with  Aa  an  incremental  change 
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in  a  (approximately  2  degrees).  Each  of  these  scans,  starting  at  A, 


result  in  a  maximum  length  Rscan(a)-  That  a  for  which  Rsean(a)  is 
maximized  is  defined  as  tne  corrected  direction  angle  for  the  line 
segment.  This  a  is  used  to  reposition  the  vidicon  assembly  on  a  new 
starting  point  and  the  procedure  repeated  until  a  new  vertex  appears 
-  coincidence  of  the  track  and  scan  no  longer  occurs. 

A  minimal  measure  of  coincidence,  or  scan  length  defined  as  a  threshold, 
serves  to  detect  the  occurrence  of  a  new  vertex  whenever  coincidence 
falls  below  thi i  threshold  value. 

The  repositioning  is  computer  as 

J'V  =  XA  +  Rmax  *cosotmax 

Xv  =  XA  +  Rmax  *sinamax 


Figure  (III -8)  illustrates  successive  stages  of  this  process  for  track 
tracing  an  event  from  point  A  to  point  A5. 


III. 9.0  THE  MAIN  PROGRAM  -  COORDINATION  OF  InH  PC/NNING  MODES 


To  demonstrate  the  data  handling  capability  of  the  main  program  digital 
simulations  of  the  analog  scanning  mode  for  the  system  were  constructed. 
Utilizing  tab  card  descriptions  of  the  tracks,  the  program  1)  Identifies 
track  entries  through  frame  edges,  2)  Traces  events  from  vertex  tc  vertex, 

3)  Assembles  event  related  data  requiring  scanning  of  more  than  one  frame, 
and  4)  Summarizes  results  of  completed  analysis  and  produces  both  card 
and  printed  copy  tabulations  of  all  vertex  coordinates  associated  with 
each  event  identified. 

The  program  utilizes  a  single  working  tape  and  instructs  an  operator  to 
mount  emulsion  frames  (card  deck  simulations)  as  called  for  by  the  logical 
assembly  of  acquired  data. 

Presented  below  is  a  description  of  the  operation  of  the  main  program. 

.he  program  will  be  described  by  following  its  operation  on  a  hypothetical 
event  which  embodies  most  possible  complexities. 

For  our  hypothetical  case  we  treat  an  event  which  covers  three  frames  of 
a  stack.  The  event  is  represented  below  by  its  projections  into  each  of 
the  three  frames  and  by  a  composite  view  of  projections  for  all  three 
frames . 

The  event,  as  sketched  below  represents  a  particle  track  entering  1  i  .  ">e  1, 
at  A.  The  track  travels  to  B,  a  vertex  of  order  2,  where  a  splitting  occurs 
forming  two  new  tracks.  One  of  these  emanating  tracks  from  the  vertex  at 
B  travels  to  D  where  it  descends  into  the  second  frame.  This  track  passes 
through  the  second  frame,  D  to  E  and  descends  to  the  third  frame  where  it 
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exits  at  edge  3,  E  to  Xj,  The  second  emanating  track  from  B  travels  to 
C  where  it  then  descends  into  the  second  frame.  In  the  second  frame  the 
track  travels  from  C  to  M,  a  vertex  of  order  2.  Two  new  tracks  are 
formed.  The  first  of  these  tracks  travels  to  F  where  it  then  descends  into 
the  third  fnme  and  exits  at  edge  2,  F  to  X3.  The  second  track  leaving 
M  travels  to  G  when  it  then  ascends  to  frame  1  and  exits  at  edge  3,  G 
to  X2. 

Flow  diagrams  detailing  the  operation  of  the  program  are  shown  in 
Gigures  (II 1-9)  through  (1 1 1 -20). 

The  main  program  begins  with  initialization  of  arrays  and  variables  which 
indicate  the  status  of  the  tracing  procedure.  As  our  hypothetical  event 
is  traced  we  will  encounter  these  arrays  and  variables  and,  therefore, 
defer  their  definitions. 

The  flow  diagram  of  Figure  (III -9),  illustrates  the  identification  of  a 
frame  (for  cur  first  frame,  FRAME=1),  the  initialization  of  pertinent  edge 
scan  arrays,  and  variables,  and  a  test  for  determining  whether  any  inter- 
frame  tracks  are  to  be  considered  in  the  frame  under  study.  At  this  point 
we  have  not  as  yet  discovered  any  interframe  transfers,  so  that  we  proceed 
to  test  for  completion  of  the  stack  edge  scan,  Figure  (III-10). 

STKFLG  and  EBGFLG  are  integer  variables  with  two  possible  values  --  1  and  0. 
S'iKFLG=l  signifies  that  the  edge  scanning  of  all  frames  in  the  stack  is 
incomplete,  while  STKFLG=0  signifies  completion.  EDGF!.G=  1  signifies  edge 
scanning  of  frame  selections  in  normal  frame  number  sequence.  EDGFLG-0 
signifies  that  the  normal  edge  scanning  of  a  frame  selection  is  interrupted 
due  to  an  excessive  number  of  interframe  track  references. 
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INITIALIZATION 


'  IGIKt.  ('ll  -9  t 


i  nit;  m  i :ai  ion  prochmikhs 


A  starting  point  for  the  edge  scan  of  the  selected  frame  is  defined  -- 
i  id ti ally  we  start  at  (0,0,0).  The  edge  scanning  routine  is  entered  and 
a  possible  tr.'ck  entry  point  is  retrieved;  namely,  point  A  on  edge  1  of 
our  hypothetical  event,  i.e. 

X0(1)  -  xA 

Y0(l)  =  Ya 

Z0(1)  -  zA 

This  entry  point  is  compared  with  the  starting  point  (0,0,G).  Xa+Y^+Za=0, 

implies  a  complete  traversal  of  the  frame  edges  and,  therefore,  a  return 
to  the  origin. 

This  ent  y  point,  then  serves  as  the  starting  point  for  subsequent  edge 
scanning,  so  that  we  redefine  XS,  YS,  and  ZS  accordingly, 

xs  =  xA 

YS  =  Ya 
ZS  =  ZA 

In  Figure  (Iil-ll),  the  entry  point  is  compared  against  previously  defined 
exit  points  --  it  may  happen  that  this  entry  point  is  actually  an  exit  point 
for  a  previously  traced  event.  Initially  the  arrays  XE,  YE,  and  ZE  are 
defined  with  zero  values  so  that  the  program  flew  proceeds  to  the  track 
entry  confirmation  procedure. 

The  distance  from  the  origin,  traveling  along  the  edges  of  the  frame,  is 
computed,  i.e. 
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EDGES  =  XA,  since  the  edge  number  equals  one. 

The  entry  point,  considered  a  vertex,  is  analyzed  to  determine  its 
multiplicity  or  the  number  of  emanating  tracks,  from  A.  The  vertex 
analysis  is  performed  by  subroutine  VERTEX.  The  arrays  ALPHA,  and 
SINBT  are  returned  to  the  main  program;  ALPHA  contains  direction  angles 
in  the  plane  of  the  frame,  SINBT  contains  the  sine  of  departure  angles 
along  the  depth,  or  z  axis,  of  the  frame.  The  first  non-zero  ALPHA (1) 
is  recorded  as  an  event,  inti  ally  EVENTS=0,  therefore,  following  Figure 
(111-12),  EVENT=1.  ORDER  is  an  array,  initialized  with  zero  values, 
which  contains  the  order  --  the  number  of  emanating  tracks  --  of  a 
designated  vertex,  i.e.  A  is  designated  vertex  i. 


Following  along.  Figure  (111-13),  the  vertex  A  is  stored  in  the  array 
COORDT.  STAR  is  the  integer  variable  denoting  vertex  numbers  of  an  eveT  t. 

With  STAR  =  1,  we  have 

COORDT (1,1)  =  X 


A  0 
COQROT(1 ,2)  =  Ya 

COORDT (1 ,3)  »  ZA 


V  IS  ENCOUNTERED 
UY  THE  TRACK  OV. 


\ 


The  order  of  the  vertex,  A,  is  determined  by  testing  thr  All  HA  array  for 
non-zero  values  which  do  not  complement  the  original  track  entering  A. 

In  general,  an  interior  vertex  is  encountered  along  a  track,  when  determin¬ 
ing  its  order  the  original  track  should  not  contribute.  For  example,  in 
the  sketch  re low . 
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DIGITAL  DATA 
STORAGE 


V  is  an  interior  vertex  --  internal  to  the  frame  boundaries  --  the  order 
of  V  is  2  not  3. 

The  order  of  vertex  A  is  1 ,  since  only  one  track  leaves  A.  The  values  of 
ALPHA(l)  and  SINBT(l)  are  stored  in  COORDA, 

C00RDA(1 ,1,1)  =  aA 

COORDA (1 ,2,1)  *  sinBA 

where  aA  and  £A  are  respectively  the  planar  direction  and  planar  departure 
angles  for  the  track  leaving  A. 

The  first  subscript  ir.  COORDA  is  the  vertex  number,  STAR,  the  second 
subscript  designates  either  planar  direction  or  departure,  and  the  third 
subscript  identifies  the  amanating  track,  presently  J=l. 

In  Figure  (III-14),  the  vertex  is  examined  for  a  point  of  exit,  if  it  is 
an  exit  point  it  is  recorded  in  the  arrays  XE,  YE,  and  ZE.  The  subscript 
S  represents  a  running  count  of  exit  points.  The  vertex.  A,  is  then 
tested  for  coincidence  --  within  a  tolerance  --  with  the  preceding  vertex. 
If  coincidence  occurs,  the  present  vertex  is  identified  as  the  preceding 
one  and  any  emanating  tracks  are  added  to  its  order.  If  coincidence  docs 
not  occur,  as  is  the  present  instance,  a  check  is  made  for  an  interframe 
transfer,  i.e.  a  penetration  into  either  an  upper  or  lower  frame.  FRAME  1 
is  an  array  which  lists  the  frames  into  which  tracks  penetrate  from  the 
event,  initially  FRAME  1  is  set  to  zero. 
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Figure  (III-14)  Test  for  Exit  and  Coincidence 


If  an  interfrair.e  transfer  does  occur,  it  is  catalogued  in  the  arrays 
FRAME  1,  COORD  1,  EVENT  1,  STARNO,  and  the  variable  FRMTRN.  In  the 
present  case,  an  interframe  transfer  has  not  yet  occurred. 

The  next  step  in  the  scanning  process  is  an  evaluation  of  the  event 
status.  Fig.  5  (III -15)  and  (1 1 1  -16)  -  The  status  of  an  event  is  evaluated 
by  comparing  the  number  of  tracks  traced  from  a  vertex  with  its  order.  If 
the  numbe1'  of  tracks  is  less  than  the  order  of  the  vertex,  this  indicates 
that  the  tracing  of  tracks  emanating  from  the  vertex  is  incomplete,  and 
therefore  the  remaining  tracks  are  traced.  The  integer  variable  START 
is  an  identifying  number  for  a  vertex  with  remaining  tracks  to  be  traced. 
The  integer  variable  COUNT  maintains  a  running  index  for  the  number  of 
tracks  emanating  from  the  vertex. 

For  the  hypothetical  event, 

START  =  1 
COUNT  =  1 

ALPHAT  =  COORDA( 1,1, COUNT)  =  aA 

SINBTT  =  COORD A (1,2, COUNT)  *  t;nSA 

X  (2)  =  X 
°  A 

Yc(2)  -  Ya 
zc(2)  -  ZA 

The  track  leaving  vertex  A,  identified  by  the  value  of  COUNT,  is  traced. 
This  track  is  traced  using  the  track  following  routine,  starting  at  vertex 
A  with  direction  angles  defined  by  ALPHAT,  and  SINBTT.  The  track  is 
followed  until  a  vertex  is  reached,  presently  vertex  B.  An  analysis  of 
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COUNT -0 
L*1 
1-1 


STAROGU)— I 


COUNT  -  COUNT  +  1 


<0  /  COUNT 
I  <\-ORDErt(ll 


START  -  I 
COUNT  -  COUNT+l 
ALPHAT  -  COOROAfl.1, COUNT! 
SINBTT  -  COORDA(l.2,COUNTI 
XO(2)  -  COOROTO.I) 
YOI2J  -  COOROTO.JI 
?0(2I  *  COORDTH.3I 


ORDER(I) 


1100  )  DATA 


Figure  (III-16)  Event  Scanning  Status 


r^‘ 


8*» 


FKilJRU  (Ill-n  -  TRACK  FOLLOWIN'. 


this  vertex,  B,  results  in  the  detection  of  three  emanating  tracks, 
defined  by: 


ALPHA(l)  =  o21 

SINBT(l)  =  sine21 

ALPHA(2)  =  a22 

SINBT (2)  =  sin322 

ALPHA (3)  =  a23 

SINBT(3)  =  sin623 

The  third  elements  of  the  arrays  represent  the  original  track,  and  will 
be  eliminated  during  the  test  for  complementation. 

Returning  to  Figure  (II 1-13),  vertex  B  is  designated  as  vertex  2,  STAR=2, 

and  its  coordinates  are  stored  in  the  array  COORDT. 

COORDT(STAR.l)  =  XD 
B 

COORDT (STAR, 2)  =  YD 

B 

COORDT (STAR, 3)  =  Zg  +  (FRAME- 1 ) *ZMAX 

The  order  of  vertex  B  is  deduced  and  stored  as  an  element  of  the  ORDER 
array  with  subscript  STAR. 

ORDER (STAR)  =  2 

The  third  emanating  track,  from  A,  is  eliminated  during  the  test  for  com¬ 
plementation  and  the  COORDA  array  is  computed. 

COORDA (STAR, 1,1)  =  ALPHA (1)  =  a2i 
COORDA  (STAR  ,2 , 1 )  =  SINBT(l)  =  sin321 
•OORDA (STAR ,1,2)  =  ALPHA(2)  =  a22 
COORD  A  (STAR  ,2, 2)  --  SINBT(2)  =  sin322 

Continuing  with  Figure  (111-14),  ■/*. rtex  B  is  tested  for  a  point  of  exit, 
which  it  isn't,  and  then  for  coincidence  with  vertex  A.  Vertices  B  and  A 


are  not  coincident,  the  program  then  checks  for  an  interframe  transfer, 
and  returns  to  evaluate  the  event  scanning  status.  Figure  (III-15). 

Since  vertex  1,A,  had  only  one  emanating  track,  which  has  already  been 
traced,  the  tracks  emanating  from  vertex  2,B,  will  now  be  traced.  The 
first  track  to  be  traced  is  defined  for  the  track  following  routine  by 
START  =  2 

COUNT  =  1 

ALPHAT  =  C00RDA(2 ,1,1)  *  a2] 

SINBTT  *  C00RDA(2 ,2,1)  =  sinB2l 
X  (2)  *  COORDT(2 ,1)  =  Xg 

Y0(2)  =  COORDT(2,2)  =  Yg 

Z0(2)  =  COORDT(2,3)  = 

The  first  track  leaving  B  is  traced  until  the  vertex  at  C  is  reached. 
This  new  vertex  is  analyzed  to  determine  its  order  and  direction  angles. 

The  coordinates  of  vertex  C,  vertex  3  (STAR=3),  are  stored  in  C.GORDT. 

COO  ROT (3,1)  =  xc 
COORDT  (3,2)  *  yc 

COORDT(3,3)  =  zc  +  (FRAME- l)*zMAX 

The  order  of  C  is  1,  since  C  is  a  vertex  due  to  the  interframe  transfer 
(Frame  1  to  Frame  2).  The  direction  angles  returned  for  tracks  leaving 
C  complement  the  entering  track  direci.l*u  angles  and  accordingly,  *  re 
deleted  from  the  list  in  COORDA. 


The  interframe  transfer,  at  vertex  C,  is  recognised,  Figure  (III-15), 
and  recorded  in  the  appropriate  arrays,  accordingly; 

FRAMEl(l)  =  FRAME+1  =  2 
C00RD1 (1 ,3)  =  0 
C00RD1 (1,1)  =  Xc 
COORD1 (1 ,2)  =  Zc 
COORD1 (1,4)  =  a21 
COORD1 (1 ,S)  =  sin621 

An  interframe  event  counter  EVENT1  is  listed  with 
EVENT1 (1 )  =  EVENT  =  1 

i.e.,  the  interframe  event  counter  links  the  transfer  to  the  event  under 
study.  The  vertex  number  of  the  event  is  recorded. 

STARNO(l)  =  STAR  =  5 

The  number  of  frame  transfers  is  updated 
FPMTRN  =  FRMTRN  +1=1 

and  the  program  returns  to  evaluate  the  scanning  status  of  the  event. 

The  evaluation  of  the  event  scanning  status  results  in  the  tracing  of  the 
second  track  leaving  vertex.  B.  This  track  is  defined,  for  the  track 
following  routine,  by 


START  =  3 


ALPHAT  »  C00RDA(2 ,1,2)  =  a22 
SINBT  =  COORDA(2,2,2)  =  sin322 
Xo(2)  =  XB 

Yo(2)  =  YB 

zo( 2)  -  ZB 

The  vertex  at  D  is  reached  by  the  track  following  routine;  it  is  then 
analyzed  and  recorded  as 
STAR  =  4 

COORDT(4,l)  =  XD 
COORDT (4,2)  =  Yd 

COORDT(4 ,3)  =  ZD  + (FRAME- l)*zMAX 

The  order  of  D  is  determined  as  1,  an  interframe  transfer  is  recognized 
(Frame  1  to  Frame  2)  and  recorded  as 
FRAME1(2)  =  FRAME+l  =  2 
COORD1 (2 , 3)  =  0 
COORD1 (2,1)  =  XD 
C00RD1 (2,2)  =  Yd 
C00RD1 (2,4)  =  a22 
C00RD1 (2,5)  =  sin622 

The  inteiframe  'ounter  is  listed  with 


EVENT1 (2)  =  1 


and  vertex  number  of  D  recorded  with 


STARNO (2)  =  STAR  =  4 

and  the  number  of  interframe  transfers  is  updated 
FRMTRN  =  ARMTRN+1  =  2 

The  program  then  returns  to  an  evaluation  of  the  event  scanning  status. 

The  evaluation  of  the  event  scanning  status,  since  all  tracks  leaving  all 
vertices  within  FRAME  1  have  been  traced,  results  in  the  tape  storage 
of  all  retrieved  data  with  an  identifying  tape  record  number.  After 
storing  these  results  on  tape,  the  arrays  and  variables 

ORDER 

COORDT 

COORDA 

are  initialized  with  zero  values. 

The  program  returns  to  the  edge  scanning  mode  and  searches  for  further 
track  entries.  For  the  hypothetical  event  under  study,  the  exit  point 
in  frame  1,  X2,  is  detected  as  an  entry  point  --  although  it  is  actually 
a  point  of  exit  with  respect  to  the  convention  established  during  the 
tracing  process.  Accordingly,  this  entry  point  is  treated  as  the  start 
of  a  new  event.  The  results  for  this  "new"  event,  i.e.,  STAR,  COORDT, 


EVENT  =  2 


the  second  event  encountered 


STAR  =  1  designates  the  vertex  X2 

COORDT(l.l)  x  component  of  X2 

COORDT(l,2)  y  component  of  X2 

COORDT(l,3)  z  component  of  X2 

The  vertex  X2  undergoes  analyses  which  results  in  the  detection  of  a 
track  leaving  X^.  This  track  is  then  followed  until  it  reaches  the  vertex 
at  G,  the  interframe  transfer  from  Frame  1  to  Frame  2.  The  evaluation  of 
the  scanning  status  for  this  "new"  event  recognizes  that  all  tracks  leaving 
all  vertices  of  Event  2  have  been  exhausted,  the  pertinent  data  for  Event 
2  is  stored  on  tape  with  an  identifying  record  number  anJ  the  program  returns 
to  search  for  new  track  entries. 

Since  there  are  no  further  entries  to  be  detected  the  program  then  selects 
Frame  2  as  the  next  frame  to  be  investigated.  Figures  (III-18)  and  (III-19). 

From  Figure  (III-18),  we  see  that 

FRAME  =  LSTFRM+1  =  2 
LSTFRM  =  2 
EDGFLG  =  1 

Since  the  new  frame  number,  2,  is  less  than  the  total  number  of  frames. 

Frame  2  will  be  mounted.  With  Frame  2  mounted  and  readied,  the  program 
first  investigates  surface  track  entries,  by  examination  of  the  interframe 
list,  and  then  searches  for  track  entries  of  new  events  originating  on 
the  edges  of  Frame  2. 


I  I  CURL  (111-19)  -  !•’ RAMI'  SHI.CCTION  IRO*-'  I  NTIiRHRAMF. 


before  continuing  with  the  trace  of  the  program  for  our  hypothetical 
event,  pertinent  information  from  Frame  1  is  summarized  in  the  table 
below. 


Event  1 


VERTEX 

NO.  VERTEX 

ORDER 

CHARACTERISTIC 

1 

A 

1 

Entry  point 

2 

B 

3 

Internal  Vertex 

3 

C 

1 

Interframe  Transfer 

4 

D 

1 

Interframe  Transfer 

Event  2 

1 

X2 

1 

Entry  Point 

2 

G 

1 

Interframe  Transfer 

Interframe  List: 

Track  Entries 

into  Frame  2 

from  Frame  1 

EVENT. 1 

EVENT  1 

EVENT  2 

FRAME! (1 )  -  2 

FRAME1  (2)  =  2 

FRAME!  (1)  1 

COORD 1 (1,1)  =  Xc 

C00RD1(2,11  =  xD 

COORDl (1,1)  = 

X^ 

U 

COORD 1 (1,2)  =  Yc 

COORD1 (2 ,2)  =  Yd 

COORD 1 (1,2)  = 

yg 

COORD 1 (1,3)  =  0 

COORD1 (2,3)  =  0 

COORDl  (1,3)  = 

0 

COORD 1 (1,4)  =  a2i 

COORD  1 (2,4)  =  a22 

COORDl (1,4)  = 

aG 

COORD 1 ( 1 ,5)=s in  21 

COORD  1  (2 ,  5)=sin322  C00RD1  (1  ,S)=sinC>G 

EVENT1 (1 )  =  1 

EVENT1 (2 )  =  1 

EVENT1 (3)  =  2 

STARNO(l)  =  3 

STARNO  (2 )  =  4 

STARNO (1)  =  1 

FRMTRN  =  2 

FRMTRN  =  1 
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In  Figure  (III-9),  Frame  2  is  identified  and  those  variables  pertinent 
to  the  edge  scanning  moae  are  initialized.  The  interframe  transfer 
list  is  examined  with  the  program  recognizing  the  existence  of  surface 
track  entries . 

The  frame  surface  track  entries  are  identified,  Figure  (III-2Q),  and 
recorded  as, 

II  =  1 

EVENT  =  EVENT! (1)  =  1 
XQ  =  COORD1 (1 1 , 1 )  =  Xc 
Y0  =  C00RD1 fl 1 ,2)  =  Yc 
Z o  =  COORD1 (11,3)  =  Zc 
ALPHAT  =  COORD1 (1 1 ,4)  =  a21 
SINBTT  =  COORD1 (1 1 ,4)  =  sim2l 

The  interframe  reference  is  then  deleted,  so  that  only  those  references 
that  remain  are  to  be  considered  during  any  further  investigations  of 
surface  track  entries,  accordingly; 

FRAME1  (II)  =  0 
FRMTRN  =  1 

The  analysis  of  the  vertex  at  C,  Frame  2,  results  in  the  detection  of  a 
track  leaving  C.  The  result  of  the  analysis  is  summarized  by 

ALPHA(l)  =  321  SINBT(l)  ■  sinful 

ALPHA (1 )  -  0,  1>1  SINBT(l)  »  0,  !>i 
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FRAME  SURFACE 
TRACK  ENTRY  IDENTIFICATION 


Y 

EVENT  =  EVENTI(II) 
START  =  STARNO'll) 
X0  =  C00RDKII.1) 
YO  -  COORDMII.2) 
ZO=  COORDKII.S) 
XO(2)  -  XO 
YO(2)  =  YO 
ZO(2)  =  ZO 

ALPHAT  =  COORDKII.4) 
SINBTT  =  COORDI(ll,5) 
FRAMEI(II)  =  0 
FRMTRN  =  FRMTRN-1 


ANALYZE  THE  VERTEX  ON 
THE  FRAMES  SURFACE 


FRAME  SURFACE 
ENTRY  PREDICTION 


EVENT  SCANNING 
STATUS 


The  non-zero  value,  ALPHA(l),  indicates  that  the  order  of  vertex  C 
is  computed  as 

ORDER (3)  =  ORDER (START)  =  2 

Continuing  in  Figure  (III-13),  the  vertex  number  of  C  is  designated 
with  a  STAR  value 

STAR  =  1 

and  its  coordinates  stored  in  the  COORDT  array 

COORDT (STAR, 1 )  =  XQ  =  Xc 
COORDT (STAR, 2)  =  Yft  =  Yc 
COORDT (STAR, 3)  =  ZQ  +  (FRAME- 1)*ZM AX 

The  tracks  leaving  C  are  tested  for  complementation,  resulting  with 

ORDER(l)  =  1 

C00RDA(1 ,1,1)  =  AT  PHA(1 )  =  a21 
COORDA  (1,2,1)  =  SINBT(l)  =  sinB2i 

C  is  tested  as  a  possible  exit  point,  which  it  isn't;  the  interframe 
test  is  passed  through  and  the  scanning  status  is  evaluated.  The 
evaluation  yields 

START  =  1 
COUNT  *  1 

ALPHA*  =  COORDA (1,1, COUNT)  =  z2\ 

SINBTT  *=  OXTOA (1,2, COUNT)  =  sinB2l 


for  the  track  leaving  C,  in  the  surface  of  Frame  2,  and  this  track  is 
traced  by  the  track  following  routine.  The  vertex  at  M  is  reached;  it 
is  analyzed  and  three  tracks  leaving  M  are  detected. 

ALPHA  (1)  =  otM1  SINBT(l)  =  sin0M1 

ALPHA(2 )  =  ayQ  SINBT(2)  =  sin6M2 

ALPHA (3)  =  aM3  SIKBT(3)  =  singM2 

The  new  vertex  is  listed,  designating  the  vertex  with  a  STAR  value,  and 
recording  its  coordinates. 

STAR  =  2 

COORDT(STAR, 1 )  =  XM 
COORDT (STAR ,2 )  =  YM 
COORDT(STAR,3)  =  ZM 

The  third  track,  defined  by  and  sinfj^,  is  eliminated  from  consider¬ 
ation  since  it  is  the  complement  for  the  original  track  entering  M. 
Accordingly , 

COORD A (STAR, 1 ,1)  =  aM1 
COORD A (STAR, 2, 1)  =  sinBM1 
COORDA ( STAR ,1,2)  =  aM2 
COORDA (STAR , 2 , 2 )  =  sintfM2 
0RDER(2)  =  3 

These  tracks  are  evaluated  by  the  program  for  their  scanning  status. 

The  evaluated  status  dictates  that  the  two  tracks  leaving  M  are  followed 


until  the  vertices,  F  and  G,  are  reached.  They  are  each  identified, 


respectively  with 

STAR  =  3,  F 
C00RDT(STAR,1)  =  Yp 
COORDT (STAR, 2)  *  Yp 
C00RDT(STAR,3)  =  Zp 


and 


STAR  =  4,  G 
COORDT (STAR, 1)  =  XQ 
COORDT(STAR,2)  =  YG 
COORDT (STAR, 3)  =  ZG 

Vertex  F  is  recognized  as  a  point  of  penetration  for  the  track  descending 
Frame  2  into  Frame  3,  and  is  subsequently  recorded  in  the  interframe 
transfer  lists, 

FRAME  1(1)  =  FkAME+1  =  3 
C00RD1 (1 ,3)  =  0 
C00RD1 (1,1)  =  Xp 
COORD1 (1,2)  =  Yp 
COORD1 (1,4)  =  aF1 
COC'RDl  (1 ,5)  =  sin8pl 

EVENT1 (1 )  =  1 
STARNO ( 1 )  =  3 
FRMTRN  =  2 
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Similarly,  vertex  G  is  recognized  as  a  penetration  point  for  the  track 
ascending  Frame  2  and  into  Frame  1,  and  is  recorded  in  the  interframe 
lists. 


FRAME1 (2)  =  FRAME-1  =  1 
C00RD1 (2 ,3)  =  ZMAX 
C00RD1 (2,1)  =  XG 
C00RD1 (2,2)  =  Yg 
COORD1 (2 ,4)  =  aG1 
COORD1  (2,5)  =  sii6G1 
EVENT1 (2)  =  1 
STARNO (i )  =  4 
FRMTRN  =  3 

After  exhausting  the  tracing  of  event  segments  from  Event  1,  that  event 
which  entered  Frame  2  through  its  surface,  the  program  then  proceeds  to 
search  for  track  entries  along  its  edges,  no  further  entries  appear. 

Following  the  edge  search  procedure  for  Frame  2,  the  third  frame  is  then 
mounted  for  investigation.  The  track  entries  into  the  surface  of  Frame  3 
vertices  F  and  G,  are  treated,  resulting  in  the  detection  of  two  tracks  - 
a  track  leaving  F  and  exiting  the  edge  at  X^  and  a  track  leaving  G  and 
exiting  at  the  edge  at  X^. 

The  program  then  searches  a^ng  the  edge  of  Frame  3  where  no  new  track 
entries  are  detected,  however  the  exit  points  X2  and  X3  were  detected 
but  only  to  De  eliminated  when  compared  against  the  list  of  exit  points. 


The  program  then  examines  all  the  data  and  recognizes  that  Events  1  and  2 
are  both  one  of  the  same,  that  is  Event  2  designated  the  track  which  left 
the  entry  point  X2  and  penetrated  into  the  second  frame.  The  program 
recognizes  that  this  entry  point  X2  is  actually  an  exit  point  for  Event  1. 

I I I. 10.0  CONCLUSIONS 

The  study  demonstrated  the  viability  of  the  conceptualization  for  a  vidicon- 
optics  hybrid  computer  system  for  a  nuclear  emulsive  track  tracer.  Major 
emphasis,  during  the  study,  was  placed  on  investigating  scanning  concepts, 
data  assemblage,  and  their  coordinated  interaction  in  a  realization  of  an 
automated  system. 

In  order  to  accomplish  the  study's  objectives  --  demonstrating  feasibility 
of  vidicon-optics  scanning  concepts  and  data  assemblage  for  particle  track 
events  --  an  optically  scaled  breadboard  model  was  constructed.  This  model, 
though  crude,  provided  the  means  with  which  to  evaluate,  experimentally,  the 
realizations  for  the  concepts  developed  during  the  course  of  this  study. 

These  experiments  verified  the  success  of  the  conceptual  scanning  regimes  -- 
edge  scanning,  vertex  analysis,  and  track  following  --  but  of  most  importance 
they  provided  insight  relevant  to  the  hardware  requirements  for  achieving 
any  useful  implementation  of  an  automated  track  scanning  system. 

It  is  necessary  to  point  out  that  the  degree  of  resolution  for  any  automated 
track  scanner  rests  almost  entirely  with  its  hardware.  Software  can  only 
treat  the  data  produced  by  hardware,  mathematical  inference  can,  of  course. 


refine  this  data;  but  the  ultimate  bound  or  resolution  is  clearly  limited 
by  the  hardware's  performance.  Our  comments  below  relate  to  criteria  for 
hardware  selection. 

The  three  major  hardware  aspects  are  the  optical  sensor,  the  optical  sub¬ 
system,  and  the  mechanical  sub-system. 

Present  day  optical  sensors,  e.g.  vidicons,  flying  spot  scanners,  etc., 
have  the  capabilities  to  offer  solutions  of  the  Older  of  500  lines  per 
inch  and  greater.  Each  of  these  sensors  should  be  considered  for  speed 
of  scanning  and  resolution. 

The  optical  subsystem  --  the  collection  of  lenses  which  project  real  images 
of  events  onto  the  face  of  an  optical  sensor  --  must  be  capable  of  provid¬ 
ing  variable  magnifications  with  minimal  distortion  and  variable  depths 
of  focus.  These  capabilities  would  provide  faster  track  following,  and 
greater  resolutions  in  determing  depth  departure  angles. 

The  mechanical  subsystem  relates  to  the  mechanisms  for  controlling  the 
motion  of  a  two  degree  of  freedom  stage  and  motion  along  the  z-axis  -- 
the  relative  motion  of  the  optical  subsystem  and  sensor  to  the  nuclear 
emulsion  under  investigation.  The  requirements  for  this  system  are 
quite  severe;  it  must  be  capable  of  resolutions,  in  its  movements,  to 
the  order  of  microns  coupled  with  fast  response  times.  These  required 
resolutions  imply  that  great  care  must  be  taken  to  mechanically  shield 
the  system  against  external  distrubrances . 


It  is  recommended  that  each  of  these  three  subsystems  be  evaluated  in 
light  of  present  day  technology  --  reviewing  existing  devices  and  systems 
including  considerations  for  their  interdependence  with  a  system. 

The  software  developed  during  this  study  is  applicable  with  minor  modi¬ 
fications  to  almost  any  realization  of  an  automatic  track  scanning  system. 
The  "vertex-to-vertex"  philosophy  remains  not  only  valid  but  appears  to  be 
a  most  efficient  scheme  for  correlating  the  vast  amounts  of  event  data. 
Most  modifications  to  the  present  software  are  envisioned  in  the  form  of 
hardware  interface  routines  and  editing  programs  which  inject  the  princi¬ 
ples  of  physics  into  evaluations  of  the  data. 


